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ABSTRACT 


A steady flow analysis is conducted for a NACA 0012 airfoil in low Reynolds 
number flows ranging from 540,000 to 1,000,000. Emphasis is placed on prediction and 
location of the separation bubble. Computational methods include the direct boundary 
layer method, the viscous-inviscid interaction method, and the time-averaged Navier- 
Stokes method. Characteristic trends in skin friction coefficient, displacement thickness, 
and boundary layer velocity profiles with increasing angle of attack are observed. 
Computational results are compared to each other and to experimental photographs 
visualizing the density flowfield using Point Diffraction Interferometry. Both the 
viscous-inviscid method and the Navier-Stokes method failed to accurately represent 
leading edge separation bubbles. The direct boundary layer method, usually considered 
of very limited usefulness due to a singularity in the underlying equations at separation, 
is shown to exhibit unexpected recovery behavior for small amounts of separation. 


Furthermore, the results near the leading edge, where separation bubbles were computed, 


were validated by the experiment. 
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I. INTRODUCTION 


While there are many reliable solution methods for high Reynolds number flows, 
the low Reynolds number regime currently has considerably fewer options. This is due 
in part to the fact that high Reynolds number flows account for most aeronautical 
applications of interest. There are, however, important applications involving low 
Reynolds numbers, such as turbomachinery blades. The more likely reason for the lack 
of reliable codes for low Reynolds numbers is the greater difficulty of accurately 
representing the flow. Most methods make use of approximations in the formulation of 
their underlying equations to obtain computational solutions in a reasonable amount of 
time. These approximations often become less and less accurate as the Reynolds number 
decreases. An obvious question which may arise concerns the value of bothering with 
such seemingly limited, simplistic codes when the state of the art is Navier-Stokes (NS) 
solvers. There are several considerations which make the effort worthwhile. First, not 
everyone has access to the supercomputers or mini-supercomputers that are necessary for 
extensive NS solutions. Even if these computers are available, user time may be limited. 
Preliminary research using a simpler method may cut down the amount of advanced 
calculations needed considerably, thus reducing overall cost. Another pertinent factor 
is the time savings. A design team for a new aircraft cannot afford months of detailed 
refinement using NS solvers when a proposal deadline looms near. Less than NS 


accuracy is certainly acceptable, especially when the significant gains in speed and cost 





reduction are considered. An efficient program only needs to calculate to the level of 
accuracy necessary to accomplish the desired goal. To this end, methods such as the 
direct boundary layer method and the viscous-inviscid interaction method offer 
opportunities to investigate low Reynolds number phenomena, such as laminar separation 
bubbles. A final point to consider is that all methods, including NS solvers, are really 
approximations in that they involve empirical models for transition and turbulence. The 
suitability of these models to low Reynolds number flows will influence the results 
directly. 

This investigation explores the various computational methods, comparing them to 
each other and to experimental data. A panel code was first developed, which provided 
some of the input needed for a direct boundary layer code. The direct boundary layer 
code was studied extensively and several modifications were made to enable further 
analysis of boundary layer profiles and transition effects. Similar calculations were then 
performed for a viscous-inviscid interaction code. Experimental interferograms, obtained 
by Point Diffraction Interferometry, provided a reliable reference for comparison. 
Finally, a Navier-Stokes code was evaluated. Emphasis throughout the investigation was 
placed on detection and location of laminar separation bubbles, as well as a thorough 


consideration of transition and turbulence models. 





i. PANEL CODES 


A. THEORY 
In potential flow theory, the flow field around an airfoil may be represented by the 
velocity potential. Considering contributions from the freestream flow and the source 


and vorticity distributions, the total potential may be constructed: 
= $+ 0,+ 0, (2.1) 
where 


o.= V.@ cosaty sina) 

o- f in ue (2.2) 

a= -f So 8 ds 
The source distributions (q) vary from panel to panel, while the vorticity strength (7) is 
assumed constant for all panels. The value of representing the flow past an airfoil by 
surface singularity distributions lies in the fact that these singularity distributions 
automatically satisfy Laplace’s equation, the governing flow equation for inviscid, 


incompressible flow: 


be: ty =0 (2.3) 


Since Laplace’s equation is a linear homogeneous second order partial differential 


equation, the superposition principle used in Equation 2.1 is valid. The boundary 





conditions include flow tangency at control points (midpoints of panels) and the Kutta 
condition at the trailing edge, requiring equal tangential velocities for the first and last 
panels. By evaluating the integrals along the airfoil surface, the potential may be 
determined at any point in the flow field. Each point is defined at a radius (r) and angle 
(9) from a chosen reference point on the airfoil. The reference point in this study is the 


leading edge. 


B. PANEL METHOD GEOMETRY 





For computational purposes, it is not feasible to evaluate every point in the flow 
field. The airfoil is represented by a number of defined points, called nodes. More 
points produce greater 
resolution and accuracy. 
One hundred to two 
hundred points are usually 
sufficient, with the larger 
numbers used for more 
complicated airfoil shapes 


or more. involved 





lations. he fi 
calculations The tines Figure 2.1 Panel Method Geometry 
connecting these nodes are 
the panels. There are (n) panels and (n+1) nodes, with the first and last node 


overlapping. Figure Figure 2.1 depicts the panel geometry. Numbering starts at the 





trailing edge, then progresses along the lower surface, leading edge, and upper surface, 
and ends at the trailing edge. The unit normal vectors (fi) are perpendicular to the 
panels and directed outward from the airfoil surface. The unit tangential vectors (t)) are 
parallel to the panels and the positive direction is defined with increasing numbering (n 
ton+1). The panels may vary in length, with the exception of the first and last panels, 


which must be equal in order to use the Kutta condition at the trailing edge. 


C. COMPUTER CODE 


1. Overview 

In order to study the steady, incompressible, inviscid flow over arbitrary 
airfoils, a panel code called panel was developed. The required input consists of the 
number of nodes on the airfoil surface, the coordinates of the nodes referenced from the 
leading edge, and the angle of attack in degrees. The program produces normalized 
velocities and pressure coefficients at each control point as output. The program was 
later modified to produce an output file compatible with bl2d, a direct boundary layer 
program described in Chapter 2. Additional input consists of Reynolds number and 
transition information. This data is not used by the program panel, but is simply 


transferred to the output file which will be used as input for bl2d. 


2. Influence Coefficients 


The use of influence coefficients leads to a straightforward procedure for 


programming the equations. An influence coefficient is defined as the velocity induced 











at a field point by a unit strength singularity distribution on one panel. For the two 


dimensional steady flow problem, the following influence coefficients are needed: 


® A‘; normal velocity component induced at the i* panel control point by unit 
source distribution on the j® panel 
nr _ 1 . 7; ij+ oe 
A a Fz sin@,-0) in — +cos(8;-6) Bl. t4] 
- "; (2.4) 
-1 i=j 
2 > 


@ A‘,: tangential velocity component induced at the i* panel control point by unit 
strength source distribution on the j® panel 


Mi 





A‘,=[sin(@,-8) 6, -cos(8,-6)In mh iej 02.5) 


=Q9 , i=j 


@ B',: normal velocity component induced at the i® panel control point by unit 
strength vorticity distribution on the j® panel 





fect ijt os ss 

B= 5—[cos(@,-0)n = -sin(8,-8)B,), i+] 2.6) 
g 

=0 ; i=j 


@ Bi; tangential velocity component induced at the i® panel control point by unit 
strength vorticity distribution on the j® panel 


B‘,= + [cos(0,-6) 8, +sin(@,-0)in-“"4], ij 
"i 


(2.7) 


2n 
1 oe 
2 > Fy 


where the geometrical quantities, depicted in Figure 2.2, are defined by: 


ry 7 y (xm,-x,° +(ym,-y 


X,+%,, 
# xm, = i “i+ 
2 
Jit diet 
ym (2.8) 
@,=arctan(2! 
je Fy 
ym:-Y;. ym,-Y, 
B, =arctan(——~2*) - arctan(—+~) 
a Aja i; ‘ 








Figure 2.2 Relationships Between Geometrical Quantities 





The formula for 8, may be verified as follows: 
a+6,+x +B; = 180° (triangle) . 
xt+b+6 = 180° (supplementary angles) 
Setting these equations equal to each other and eliminating common terms, 
Bi; = b-a 
Inspection of the diagram shows that angle b is in fact the arctangent of the quantity in 
parentheses in the first term of the formula for 8. Likewise, angle a matches the second 


term. 
3. Program Description 
a. Boundary conditions 
The first boundary condition requires flow tangency at control points: 
(V"%,=0 , i=1,2,...7% (2.9) 


In terms of influence coefficients (with V,.=1), 


yet A" aly yr B® ,+sin(a-8)=0 , i=1,2,....0 (2.10) 





The second boundary condition is the Kutta condition, which states that 


the pressures on the lower and upper panels at the trailing edge must be equal if the flow 


is to leave the trailing edge smoothly. Using a form of Bernoulli’s equation, 








oF V, 
C,= PrPe 2 -(—alyp (2.11) 
Yep ve Vv, 
the pressure equilibrium also implies equal velocities for incompressible flow. Since the 


normal velocities are taken to be zero, the boundary condition may now be stated as: 


(V),=-V, (2.12) 
where the negative sign is strictly due to the adopted convention of positive tangential 
velocities in the direction of increasing node numbering. Since the flow is positive to 
the right (as shown in Figure 2.1), the panels downstream of the front stagnation point 
will have negative values for computational purposes only. It is important to note that 
not all the lower surface panels have a reversed sign, only those downstream from the 
stagnation point. This is especially significant for non-symmetrical airfoils or any airfoil 
at an angle of attack. 


In terms of influence coefficients, the normalized equation becomes: 


-YY., A'yad-v 4, BYy-cos(a-6,) = e135 
Dy, H'yq) 7D, Bi, +008(a-8,) 
b. Solution procedure 
Equations 2.10 and 2.13 represent a linear algebraic system of (n+1) 


equations and (n+ 1) unknowns. The unknowns are the source strengths which vary from 


panel to panel (q,...q,) and the vorticity strength +. 











Expanding and rearranging Equation 2.10 for an example airfoil of 


n=73 nodes and panels results in: 


AMG + AMG te tV(B%, + By t+ BaF -sin(« -8,) 
A Gd, + AnnG tt YB "ay + By te. B", 43) = -sin(a-6,) (2.14) 


A' 191 tA tt YB "a1 + Bigot tB "73,73) = -sin(« -8,,) 


The equations now readily lend themselves to solution in matrix form. Recasting with 
a simpler notation, the A*, terms (coefficients of q,) may be renamed a; and the sum of 
all B°; terms in parentheses (coefficients of y) renamed a,,,,, where i=1,2,...,n and 
j=1,2,...,n. The terms on the right sides of the equations may be renamed b,. 

The (n+ 1)* equation, or in this example, the 74° equation, comes from 


Equation 2.13 in a similar manner: 


A‘ +4 9g )G +A‘ tA "Got (A's 754A 99.73) 473 + 
Y((B 4, +B og) + (Bot B32) + ~ +B 73° B aaa3)] = (2.15) 
cas (a -6,) -cos(a -6,,) 
The coefficients of q; may be renamed a,,;. All of the B' terms in the brackets together 


form the coefficient of y, now renamed aj474,. The entire right side of the equation 


constitutes the new term b,,. 


10 


Finally expressing this system in a concise matrix form for the general 


case, 


By By yg were 4.1 | (9, b, 
4, 9, G3 + Bonet | |Fo b, 
a a ore iC ae b, 2.16) 
G1 Any gg ssrenenene Ont 11a b, 
Batt Fn Fneig or Fnatpot| LY bet 


This system is solved in the program using a Gaussian Elimination subroutine. 
With the values of the q; and y known, the velocity at each panel 


control point may be calculated: 


V, =o [A‘q]+cos(a-6) , i=1,2,....n (2.17) 
The total velocity is equal to the tangential velocity due to taking the normal velocity to 


be zero. 


c. Numerical techniques 
Although programming most of the described procedure is a relatively 
straightforward task, there are a number of potential pitfalls. Different programming 
languages each have their own special rules and format, but the following hints for 
FORTRAN, a commonly used engineering language, apply to many others as well. 


All angles entered into the equations, such as a, 8, and 6, must have 


values in radians. The angle of attack (a) is an input parameter that should also be 











converted relative to the zero lift line of the airfoil in order for the resulting velocities 
to. match published values. The angle @ is used in calculating the influence coefficients. 
@ is the angle of a panel from the positive x axis, counter-clockwise positive. The 
formula given for 6, involving taking an arctangent, will produce the correct physical 
angle if interpreted correctly. Many programming languages use -2/2 to 2/2 as the 
default range for the standard inverse tangent function, which uses only one value for an 
argument. The function cannot determine whether a negative sign was in the numerator 
or denominator. A problem arises when an angle is actually in the 2°? or 3° quadrant 
because the function will assign values from the 1" or 4° quadrant. The values for @ 
must be in the range -7<6<- to work properly. This may be accomplished by using 
the ATAN?2 function in FORTRAN, which accepts both a numerator and denominator 
as arguments and assigns quadrants correctly. 

The angle 8 may be calculated from two inverse tangents, as presented 
in the formula. However, a more efficient algorithm may be used for computer 
calculation. Recalling that 6 = b-a, it follows that: 

tan 8 = tan(b-a) 


= sin(b-a) 
cos(b-a) 


= sin b cos a - cos b sina 
cos b cos a + sin b sina (2.18) 


The sines and cosines for the angles a and b may be easily determined from the geometry 
of Figure 2.2. For example, sin b = (ym; - y;,,)/hypotenuse b. Noting that all the 


hypotenuse values may be cancelled out in Equation 2.18, the arctangent of the angle 


12 








may be found with simply x and y differences. The ATAN2 function in FORTRAN 


i 


should also be used to calculate 8. 


D. RESULTS 
1. Eppler E585 Airfoil 
The first investigation was conducted for an Eppler E585 airfoil (shown in 
Figure 2.3), with n=71 nodes (or panels). It is one of a series of airfoils without flaps 
designed for the Reynolds number range of sailplanes, about 100,000 to 500,000. The 
angle relative to the zero lift line is 
5.53°. If not compensated for, the results 
would all be shifted by this amount. 
Figure 2.4 shows the results of the 
program panel for « = 3, 7, 11°. The 
published data in Figure 2.5 [Ref. 1] 
compares favorably. The velocities match 


the given distribution well except for 





slight deviations at the trailing edge. This Figure 2.3 Eppler E585 Airfoil 
difference can be attributed to the higher 

order panel method used in Eppler’s computations. The panels are defined by third 
degree polynomials whose coefficients are computed by a spline program. Additional 


points are splined in as needed. For the best precision from a panel method, the steps 


must be smaller near the leading and trailing edges. The simpler method of connecting 














« relative to the zero-iift line 





Figure 2.5 Velocity Distribution Computed by Eppler 





just the given airfoil coordinates naturally leads to less accurate results in the critical 
trailing edge area. Experimentation with linearly interpolated extra points verified that 
a better resolution could be obtained in this manner. For most cases, the straight line 
segment panels using only the given coordinates produce quite sufficient resolution. with 
a great advantage in computational speed. 
2. NACA 0012 Airfoil 

This is a well known airfoil (shown in Figure 2.6) used for many studies and 
comparisons, elementary to advanced, as well as practical applications. It is one of the 
original 4-digit series of 1932, 
where the first two digits indicate 
camber amount and location (00 is 
symmetrical) and the last two 
digits indicate maximum thickness 
in percent mean aerodynamic 
chord (12%). Figures 2.7 and 2.8 
show a comparison of program Figure 2.6 NACA 0012 Airfoil 
results and those provided by Anderson [Ref. 2] for the NACA 0012 airfoil at 9° angle 
of attack. The pressure distributions both come to a suction peak of -5.2. It should be 


noted that Figure 2.8 differs slightly from the original reference plot in that the error in 


the decimal place of the ordinate values has been corrected. 
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Figure 2.8 Pressure Coefficient Provided by Anderson 
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E. USER’S GUIDE 

These detailed instructions are accurate at the time of writing, oriented toward the 
Advanced Computation Laboratory (rm. 136) of the Naval Postgraduate School. 
However, due to the dynamic nature of any computer environment, some instructions 
may change over a period of time. If difficulty is encountered, check with the manager 


for changes in: 


i. the account location of the files 
2. plotting package availability or procedures 


3. computer informal addresses (i.e. madmax, suzqt, indigol, etc.) 


In all instruction lines, the prompt which appears on the screen is in normal type, while 
literal user input is in boldface. User input which is general and requires the actual 
word instead will be in isalics. For example, filename could be replaced by vel.dat. 
Although not specifically stated on each line, a carriage return (ENTER) is assumed. 


As many commands are case sensitive, type each line in the format shown. 
1. Stardent 
a, Program operation 
Using the Stardent terminal, login and change to the directory to be 


used for the panel code, making a new one if desired (i.e. mkdir paneldir, cd 


paneldir). Then type: 


> cp /alpha/acctname/panel panel 








to copy the compiled program from the account where it is stored. Also copy a sample 
input file fora NACA 0012 airfoil: 

> ep /alpha/acctname/points.dat points.dat 
These two files are all that is needed to begin. To run the program, type: 

> panel 
Answer the prompted questions, using 100 points and 0 for the transition code. When 
finished, the program will respond with: 

Calculations complete, output files are: 

vel.dat, cp.dat, bl2d.dat 
The first two files are simply x-y type column data for plotting the velocity and pressure 
distributions. The last file is created for use with another program, bi2d, described in 


the next chapter. 


b. Plotting procedures 
To view the resulting velocity distribution, the file vel.dat can be used 

with any standard plotting package. On the Stardent, gnuplot can be used. First type: 

>xterm -t & 
to open up a tektronics window for plotting that will operate in background. When it 
appears, move the cursor there and type: 

> gnuplot 
Note that it is somewhat "messy" to work in this window because mistakes do not 
disappear from the screen with a backspace. The computer will take the overwritten 


characters as the input, however. Control-C may be pressed instead to simply type the 
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line over or type clear to erase the whole screen when needed. At the gnuplot prompt, 
type 

gnuplot > set term tek40xx 
to set the proper terminal type for plotting. To plot the x-y format data file vel.dat, 
type: 

gnuplot > plot "vel.dat" 
For options in gnuplot such as scaling axes and plotting multiple data sets, see Appendix 


B. 
2. IRIS 


a. Using a Stardent window 
To simply operate the program panel from an IRIS terminal, login and 
open up a Stardent window using the following procedures. First, click the right mouse 
button inside the original console window and select Clone. Move the cursor to the new 
window and change to the desired directory on the IRIS, making a new one (mkdir) if 
desired. Then type: 
> telnet suzqt 
Login to the Stardent and follow the instructions in the previous section to operate the 
program. Do not use the same plotting procedures when complete, however. Keeping 
both windows open, move the cursor back to the IRIS window and type: 


>rep suzqt:/alpha/loginname/directory/vel.dat vel.dat 
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This command remote copies from one system to the other. The loginname on the 
command line is the name of the account, usually the user’s last name. The directory 
is the one created on the Stardent where the program was run. This command will only 
work if an account is held on both the Stardent and IRIS under the same loginname. If 
this is not the case, files can be transferred using the file transfer utility ftp, described 
in Appendix B. 
b. Using the IRIS 
It is also possible to do all calculations directly on the IRIS. This may 
be more useful when a user holds an account on the IRIS but not on the Stardent. After 
obtaining a copy of the source code panel.f from the Stardent using the ftp procedures 
(Appendix B), compile it for the IRIS: 
> £77 -O3 -o panel panel.f 


Program operation is as described for the Stardent. 


c. — Plotting procedures 
The plotting package usually used on the IRIS is XYPLOT. At the 
prompt, from the directory with the plotting data (such as vel.dat) type: 
> xyplot 
Answer the questions that follow: 
Name of Ist input file? vel.dat 
Name of 2nd input file? (press ENTER since only 1 file to plot) 


Default configuration file? | (ENTER, since none specified yet) 


20 





A new window will pop up with the plot. The mouse can be used in this window to 
change the appearance of the plot in many ways. The plotting program is very user- 
friendly and can be operated with little prior instruction. 

Gnuuplot is also available on the NPS IRIS. It is not necessary to open 


a special window for the plot because one will be created automatically when the 


program is invoked. 





Il. DIRECT BOUNDARY LAYER CODE 


A. THEORY 
The two-dimensional flow around an airfoil may be represented by the simultaneous 


solution of the continuity equation, 


9p +5Pu , OPV_p G.1) 
a ax oy 


the components of the momentum equation, 


Spu 2pu? , Spuy --P 20m on +pf. (3.2) 
a ae Cay ay 
Sov , dour , dev? __& , xy , Oy og G.3) 
em yl ey 


and the energy equation, often collectively referred to as the Navier-Stokes equations, 
although technically this name applies only to the momentum equation applied to a 
Newtonian fluid. Equations 3.1-3.3 are expressed in general form for unsteady, 
compressible, viscous flows with body forces. In practice, however, such a complete 
solution is not usually feasible, or even possible, for many applications. The 
development of the Thin Shear Layer (TSL) equations, sometimes referred to as the 
boundary layer equations, enables a computationally practical scheme for solving the 


flow. 
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Considering the steady, incompressible, flow around an airfoil with no body forces, 
some initial simplifications may be made. The energy equation becomes decoupled from 
the other equations and is no longer needed in the solution. The continuity equation 


immediately reduces to 


ou , Ov (3.4) 


—+—=0 


ax oy 
Eliminating the time dependent term, expanding the partial derivatives, subtracting u 
times the continuity equation, and dividing through by p in the x momentum equation 


yields 
Ai yO. 1 3,1 Wm, 1 ry G5) 
with a similar expression for the y component. The equation is further developed by 


using the assumption of a Newtonian fluid, in which stress is proportional to rate of 


strain: 


Ou 

=2n( 4 3. 

ag 2 a (3.6) 
ou, Ov 

Oy Mata? (3.7) 


where yp is the viscosity. In a constant-property flow, y may be taken outside the 
derivative after substituting Equations 3.6 and 3.7 into Equation 3.5, and may be 


rewritten in terms of the kinematic viscosity, »=y/p. Thus, the x-component momentum 
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equation for a Newtonian fluid in constant-property, steady, two-dimensional, laminar 


flow is 


4.4 yh. 1 4 Fu, Puy (3.8) 
a ay pa a? ay? : 


To incorporate the effects of a turbulent flow, all instantaneous flow quantities are 
replaced by a mean term plus a fluctuating part, i.e. u=utu’. Expanding, eliminating 
zero products, and rearranging the equation gives 

2 TT 
= Fu, Pu _ou'? du’ (3.9) 
a ay pa ax? ay* a oy 


where the overbars on the mean velocity components u and v are omitted for simplicity. 





The extra turbulent stress terms in Equation 3.9, as compared to Equation 3.8, are often 
called the Reynolds stresses. 

To obtain the TSL equations, an order of magnitude analysis is applied with the 
assumption of 5/1<1. In other words, the boundary layer (of thickness 6) is very small 
compared to the characteristic length of the body. Using the following order of 
magnitude approximations 

Ba he he (3.10) 
7 do 8 ax i 
where the subscript e refers to the edge of the boundary layer, the first term in the 
parentheses and the first Reynolds stress in Equation 3.9 may be neglected. When 


compared term by term to the x-component equation with the assumption u>v, the y- 
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component of the momentum equation reduces to the approximation of constant pressure 
in the normal direction. Summarizing, the two-dimensional, incompressible, steady 


boundary-layer equations for both laminar and turbulent flows are: 





ou wv 
Gh Wg (3.11) 
dy 
ye ey 1B, y Fh UY (3.12) 
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a Gy pax dy? ay 


oP . (3.13) 


ay 
Note that these equations are for a surface coincident with the x-axis. In order to use 
the equations for an airfoil, the airfoil surface must be "unwrapped" onto the x-axis. The 
usual x/¢ and y/c coordinates which define the airfoil must be transformed to a surface 
coordinate. 

The applicable boundary conditions on the surface of a solid airfoil are 
y=0 u=0 v=0 

and at the outer edge of the boundary layer, 


y=5 — u=u.(x) 
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B. COMPUTER CODE: BL2D 


1. Overview 

The program bl2d, developed by T. Cebeci, provides a solution to the 
boundary layer equations. The same restrictions apply to the program as to the equations 
on which it is based: two-dimensional, steady, incompressible, viscous flow. The 
program accepts input of Reynolds number and prescribed transition locations, as well 
as panel coordinate and velocity information computed by a separate routine. Output is 
generated for many features of the resulting boundary layer, including skin friction 
coefficient and displacement thickness. Run time is less than one minute on a Stardent 


computer and less than two minutes on a personal computer (PC). 
2. Models 
a. Turbulence Model 
In order to use equation 3.12, an expression must be found for the 


Reynolds shear stress term. Since it is not feasible to attempt calculating the actual 


value, empirical models are usually used. One such model is the eddy-viscosity concept: 


“8H =e Se (3.14) 


where €,, is an empirical term called the turbulent eddy viscosity. Another model is the 


mixing-length concept first proposed by Prandtl in 1925: 
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where | is the empirically determined mixing length. Although both of these models 
i were originally derived based on erroneous physical arguments, they have nevertheless 
produced remarkably successful results for many applications. A third model, which 
incorporates the eddy-viscosity concept, is the Cebeci-Smith (CS) model, in which the 
viscous region is divided into an inner layer and an outer layer, each with its own 


formula for e,: 


Em | _ Oy 2 
= SO leyRe lice Yn’vy,, 


(3.16) 
[s}- =0.0168 /Re, In, fly, 
where 
u 
wea 
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(3.17) 
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G= i Me e 2Re -134 
1200\u_} * 
and the Falkner-Skan variables , n, and f are used. The term +, is a factor which 


models the length of the transition region, explained further in the next section. This 
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turbulence model is used in the program. Using Equation 3.14 allows rewriting the 
momentum equation for a turbulent flow, Equation 3.12, in the same form as a laminar 


flow: 


14h gy St. 1D , O (4 Oy (3.18) 
& dy pa ay ay 


where b=y+e,. Thus, the only computational difference between laminar and turbulent 


boundary layers is the addition of the turbulent eddy viscosity. 


b. =‘ Transition Model 

The determination of the location of transition from laminar to turbulent 
flow is one of the most critical factors in the success of many computational efforts to 
predict or reproduce physical phenomena. Yet, even in today’s age of supercomputers, 
this area of research remains widely neglected. The most advanced Navier-Stokes 
solvers often ignore the issue entirely, arbitrarily declaring the entire flow to be 
turbulent. Others make the effort of computing a transition point, at which the flow 
instantaneously changes from laminar to turbulent. These may be quite reasonable 
approximations for many applications, especially when the Reynolds number is high. 
However, there are still a number of important flows that require more accuracy. Until 
the transition mechanism and the many, varied factors which can affect it are more fully 
understood, programmers must rely on the traditional engineering approach of modeling. 

The program incorporates a transition model determined by Chen and 
Thyson, utilizing a transition range rather than a point. This range is neither fully 


laminar nor fully turbulent. It is a region of intermittency, in which turbulent spots 
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gradually appear with progression in the streamwise direction. This feature has. been 
shown to be essential for low Reynolds number flows. The convergence of the entire 
boundary layer solution is very sensitive to transition related factors, such as the input 
value supplied for the start of transition. An experimental value may not have been 
measured accurately and an empirically calculated value may deviate from the actual 
onset of transition. If the code does not run and all other input has been verified to be 
correct, a solution may often be obtained by experimenting with slight deviations in the 
transition location for the upper surface specified as input. The lower surface value does 
not exhibit the same sensitivity. 
3. Numerical Techniques 
a. Transformation of Airfoil Coordinates 

The x/c and y/c coordinates of the airfoil are supplied as part of the 

input, called xc and yc in the program. Starting from the stagnation point, the program 


redefines these coordinates into a single parameter corresponding to a surface distance: 


HRA-aty (xe,-x0,_.) +(9¢,-Ye,_4)° ed 


Thus, the variable x used internally by the program in the boundary layer equations is 
this surface coordinate. The values are printed under the heading (S) in the output. 


b. ‘Transformation of Variables 


The well known Falkner-Skan transformation is used to transform the 


variable y: 








ue (3.20) 


where y is the normal coordinate along which the thickness of the boundary layer is 
measured. The dimensionless similarity variable 7 eliminates the growth of the boundary 
layer in laminar flow and reduces it in turbulent flow. This enables larger steps in the 
streamwise direction and improves computational efficiency. The x transformation is 
simply a scaling by the reference length, usually the chord for a airfoil, so that £=x/c. 
Since the surface distance x is used, these will not be the same as the input x/c 
coordinates. 
The dimensionless stream function f(x,y) is defined by 


Y 


fu,vx 





Sm) = (3.21) 


Equations 3.11 and 3.18 and the boundary conditions may be rewritten in terms of the 


new variables: 
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where », is the transformed boundary layer thickness corresponding to 6 and R, is the 


Reynolds number based on reference velocity u.. and reference length, the chord for an 








airfoil. The prime denotes differentiation with respect to 7. The dimensionless pressure 


gradient parameter m is defined by 


(3.23) 








The velocity components v and v are related to the dimensionless stream function by 


f_d F 

u=t va-fuvx = fast 3.24 
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c. Keller Box Method 

Equation 3.22, a second-order partial differential equation, may be 

solved by various numerical methods such as the Crank-Nicholson or Keller Box 

methods. The latter method, depicted in Figure 3.1, has proven to be efficient for 


boundary layer calculations. The Keller Box method first requires reformulating higher 





Figure 3.1 Grid Box for Centered-Difference Approximations 
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order equations into a set of first order equations. At each rectangular grid section, these 
equations are approximated using centered-difference derivatives, averaging values at the 
center of the "box". The truncation error is second order. The resulting implicit, 
nonlinear difference equations are linearized by Newton’s method and solved by a block 
elimination method. 


Using the following definitions, 


fizu u’=v (3.25) 


Equation 3.22 may be expressed as a first order system: 


+ (n+) + ~ gf 
(bvy' ty m(1-u?) = bus "8 (3.26) 


n=0 u=0 f=fQ) , nen, ul 
Note that the u and v in Equations 3.25 and 3.26 are not the velocity components. They 
are two new, arbitrarily selected names of variables for the expression of the first order 


system. Since the program was coded using these particular variable names, the present 





numerical discussion will use them for consistency. 


Referring to Figure 3.1, the grid points may be described as 


graertsk* on=1,2..,N  &=0 
nyermyrth,  J=1, 20S Hy =0 (3.27) 
nsFN. 


where the capital letters N and J are the maximum number of grid points used in the 
streamwise and normal directions, respectively. The superscript n is not an exponent, 


but a counter just like the subscript j. This upper and lower notation allows using both 
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counters on the same variable when needed. Considering one box of the grid, the finite 
difference approximations of Equation 3.25 may be written for midpoint of the right side, 


segment P,P,, using centered difference derivatives: 


ie i “Fi up tuys 


ft ay", 
Bi = ; (3.28) 
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h 2 ie 


Equation 3.26 may be approximated in the same manner for the midpoint of the box by 


centering first in one direction then the other. The resulting finite difference equation 
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The boundary conditions of Equation 3.26 are rewritten at £=£° as 


fof, u=O0  up=l (3.30) 
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d. Newton’s Method 
Equations 3.28, 3.29, and 3.30 comprise a set of 3J+3 equations and 
33-+3 unknowns (f?, u, v, where j= 0, 1, 2, ..., J), with ff, uf", and v>? known. 
Newton’s method is applied to linearize this system. The method assumes that an 
approximate solution is known, either from the preceding iteration cycle or from the 


previous streamwise station. Then small unknown quantities are added to the 





approximate solution. Using the arbitrary iteration variable i (the superscript n omitted 


for clarity), 


fi =fis af! ui =u bu; yi =v} av} (3.31) 


with i=0 corresponding to known values at the previous streamwise station (£*"), these 
expressions may be substituted into Equations 3.28 and 3.29 for the unknowns. After 


dropping higher order terms of 4, a linear system of equations results: 
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where the right hand sides are 
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and the coefficients are 
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The boundary conditions of Equation 30 become 


Bfp=0 bu-0 bu; =0 (3.33) 
These equations may be easily identified in the subroutine COEF of bl2d. Since they 
may be arranged into a block tridiagonal structure in matrix-vector form, the subroutine 
SOLVE uses the efficient block elimination method to solve for the small 6 quantities. 
The iteration of Newton’s method continues until the small quantities are small enough 


to be neglected. 


4. Program Modification for Boundary Layer Profiles 
At each station along the airfoil surface, the program calculates the u 
velocities for each value of 7 in the grid. Eta is the coordinate in the normal direction 
representing the transformed boundary layer thickness. In order to retrieve the physical 


boundary layer thickness, an inverse transformation is required: 
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7 (3.34) 


Since the kinematic viscosity shows up only indirectly in the non-dimensional 


form of the Reynolds number, the actual equation used is 


y=n | — 3.35) 





where R, =pu,x/p, or equivalently, u,x/p, and x and u, are used in the non-dimensional 
forms of x/c and u,/u,,. The value of x used here is the surface distance. 

Plotting the shape of the actual velocity profile at a given station requires the 
station number (NX), the corresponding x coordinate, the u velocities, and the 
corresponding y values. All of these values are provided by the original program or 
Equation 3.35. Velocity profiles may be plotted at this point, but only the shape will be 
revealed. To visualize the growth of the boundary layer, the height of the boundary 
layer at each station is needed. This may be determined by finding where the u velocity 
has reached freestream velocity, indicating the. edge of the boundary layer. 
Computationally, this is accomplished by allowing u to reach 0.995 of u,., the freestream 
velocity. Even though the remaining u values in the grid above this height will still be 


calculated by the program, no more values are written to the plotting output file. 
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5. Program Modification for Estimating Transition Location 
The original program uses input values to specify the onset of transition. The 
transitional flow region is then calculated using the Chen-Thyson model, shown as ‘yj, in 
Equation 3.17. In order to provide an initial estimate for the transition location when 
no other method of determination is available, a modification using Michel’s criterion 
was incorporated: 


Ry =1.174(1 +20) Re, 8 (3.36) 
Xp 


where R,, =u,0/v is the Reynolds number based on momentum thickness at transition and 


Re,, is the Reynolds number based on the transition location. 


C. RESULTS 


Studies were conducted to: 


@ validate the program by comparing to known data 
® determine the effect of prescribed onset of transition 


@ investigate the possible occurrence of zero or negative skin friction to indicate the 
start of a laminar separation bubble before breakdown of the code 


® analyze boundary layer velocity profiles on the airfoil upper surface 
@ evaluate the effects of changes in Reynolds number 


® assess the validity of obtaining an unsteady boundary layer solution by extracting 
steady velocities from unsteady pressure distributions 
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All studies presented are fora NACA 0012 airfoil defined by 100 points. The Reynolds 
number is 540,000, except where noted in the validation study and the Reynolds number 


effect study. 


1. Program Validation 

In order to ascertain that the results of the program could be considered 
reliable, an initial test case was run to compare with published data [Ref. 3]. The test 
conditions were an AOA of 0° at a Reynolds number of 6,000,000. The first results 
obtained were in the expected range but the curves were not smooth. The problem was 
traced to a very smal! discrepancy in the original airfoil coordinate input file provided 
with the program which was not noticeable when the file was checked by plotting. After 
generating new airfoil coordinates and running them through panel to get new velocities, 
smooth boundary layer results were obtained. Figure 3.2 shows the computed skin 
friction coefficient and displacement thickness and Figure 3.3 shows the published 


results. The plots exhibit excellent agreement. 


2. Transition Onset Location 
In the next study, the convergence of bi2d at a lower Reynolds number of 
540,000 was investigated. The input value for the location of the start of the transition 
range was found to be the most critical factor. Starting at an angle of attack of 0°, 
transition values obtained from incompbl (Chapter 4) were initially used. The program 


converged to a reasonable solution as determined by viewing plots of the various output 


files. At 2° however, the program would not converge using the estimated transition 
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Figure 3.3. Reference C, and 5*, NACA 0012, AOA=0°, Re=6,000,000 
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value. Experimentation with this input parameter showed that moving it back usually 
made the convergence even worse, but in moving it forward, a point was reached where 
a solution could be obtained. Moreover, this solution was one that could be reasonably 
expected based on the previous solution at 0°. Similar experimentation was performed 
in AOA increments of 2° until excessive separation at high angle of attack caused the 
code to break down. The same study was also conducted later using the version of bi2d 
modified to make transition estimates. Table 3.1 shows a comparison between the values 
estimated by the two programs and the value actually needed for the program to 


converge. 


Table 3.1 COMPARISON OF TRANSITION ONSET 





The first significant point to note is that the transition values produced by the 


modified version of bi2d are consistently close to the values produced by incompbl, 





showing that the criterion has, in fact,.been programmed correctly. The differences can 


be attributed to several factors: 


@ slightly different values are input into the criterion equation for each program 


© incompb! outputs the x/c value of the nodal point nearest the calculated transition 
onset, as opposed to the actual value 


@ if the onset of transition is located inside a separation bubble by the initial 
calculation, incompbl arbitrarily moves it to the start of the bubble 

The third column shows how far forward the transition point was moved to obtain 

convergence, where the first AOA did not exhibit sensitivity. The last column shows the 

percent difference between the best estimate and the value required for convergence. 

Most were fairly close, with even the worst case less than 30% forward of the first 

estimate. This sets a reasonable bound for necessary experimentation with the transition 


location. 
3. Laminar Separation 


a. Skin Friction Coefficient and Displacement Thickness 
Figures 3.4-3.10 show the progression of skin friction coefficient (C,) 
and displacement thickness (6*) as the AOA ranges from 0 to 12 degrees. The transition 
onset may be observed as the point where C, reaches as minimum then dramatically 
increases, indicating the change from laminar to turbulent flow. : The transition point 


moves forward as the angle of attack increases. The minimum value of C, decreases 


with increasing angle of attack. When the C, reaches zero, separation is indicated. 








Figure 3.5 BL2D: C, and 6+, NACA 0012, AOA=2°, Re=540,000 
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« Figure 3.7 BL2D: C, and 6*, NACA 0012, AOA=6°, Re=540,000 
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Figure 3.9 BL2D: C, and 6+, NACA 0012, AOA=10°, Re=540,000 
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Figure 3.10 BL2D: C, and 6*, NACA 0012, AOA=12°, Re=540,000 

Knowing that the boundary layer equations break down when separation 
occurs, the anticipated information of this study was the x/c location of transition for as 
many angles of attack as possible before the C; became negative, where it was assumed 
the program would not run. A beneficial discovery of this study is the ability of the 
program to recover from mild amounts of separation with meaningful results. At 4°, the 
first encounter with separation may be observed as the C; just dips below zero. The 
pattern over the remaining airfoil surface suggests a separation “bubble” after which the 
flow reattaches, as opposed to near-stall separation. The program exhibited this recovery 
behavior all the way to 10°. The final plot at 12° shows that even though a solution was 
produced, convergence was not attained and the results were meaningless due to the 


greater amount of separation. 





b. Boundary Layer Velocity Profiles 

To complement the skin friction study and to further investigate the 
pattern of laminar separation, the program was modified to calculate and produce plotting . 
output for velocity profiles at evenly spaced intervals along the top surface of the airfoil. 
Results are shown in Figures 3.11-3.16 for an AOA range of 0 to 10 degrees. Boundary 
layer growth is evident as the angle of attack increases. Furthermore, the region most 
prone to separation, as revealed by the point of inflection in the velocity profile, moves 
forward with increasing angle of attack, confirming the indications of the skin friction 


plots. 
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Figure 3.11 BL2D: Velocity Profiles, NACA 0012, AOA=0°, Re=540,000 


46 







Station 0 is 


stagnation point 


6 15 20 25 30 35 40 45 
Airfoil Upper Surface Station 





Station 0 is 


stagnation point 


[oan EE > EE o> EE <> <> > > 


“6 15 20 25 30 35 40 45 


Airfoil Upper Surface Station 
Figure 3.13 BL2D: Velocity Profiles, NACA 0012, AOA=4°, Re=540,000 
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Figure 3.15 BL2D: Velocity Profiles, NACA 0012, AOA=8°, Re=540,000 
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Figure 3.16 BL2D: Velocity Profiles, NACA 0012, AOA=10°, Re=540,000 





4. Reynolds Number Changes 
The direct boundary layer code was run for Reynolds numbers of 540,000, 
750,000, and 1,000,000. A representative sample of the results at 2° is presented in 
Figure 3.17, showing a comparison of the skin friction coefficients. As the Reynolds 
number increases, two effects may be observed. First, the transition point moves 
forward, so there is more turbulent flow. This occurs on both surfaces, although only 


the upper surface is plotted. In addition, the minimum value of C,; increases. Both 





effects suggest that separation is less likely to occur as Reynolds number increases, all 


other conditions being constant. 











u000000T=88u 


unQ000SL=8%u 
«0000bS=2uu 









uz uoTtIeTazAeA “Hep Z=vVOV Z100 YOWN 











Effects of Reynolds Number Changes 





Figure 3.17 BL2D 


5. Unsteady Boundary Layers 

By correlating an unsteady pressure distribution with a similar steady one, the 
question arises as to whether the steady direct boundary layer code can process the data 
in the form of a steady velocity distribution and produce results that correspond to the 
unsteady case. In the studies of Neace {Ref. 4], it was shown that an unsteady pressure 
distribution at a certain angle of attack has a closely matching steady pressure distribution 
at an angle of attack approximately 2° lower. For example, 10.14° unsteady corresponds 
to 8° steady. An unsteady panel method program called U2DIIF was used to calculate 
unsteady pressure distributions for ramp motion. With the hypothesis that the boundary 
layer characteristics, steady or unsteady, are driven by the pressure distribution, the 
pressures were converted to "steady" velocities for the unsteady angles of attack using 
Bernoulli’s equation. Steady velocities for steady angles of attack were computed with 
the program panel. In Figure 3.18, the comparison between steady and unsteady 
velocity distributions reveals excellent agreement on the upper surface and fair agreement 
on the lower surface. Since velocity is the main input to bl2d, and the upper surface is 
much more critical (for positive angles of attack), the strong correlation suggests that the 
steady boundary layer code may indeed be able to produce results for the unsteady case. 
Many attempts were made to obtain such results using all of the methods explained in 
other sections to facilitate convergence. No solution was found for any of the cases 
investigated, however. Apparently, the differences on the lower surface had more effect 
than anticipated. Future investigations could include a modified velocity input, where 


the unsteady: upper surface is spliced with the corresponding steady lower surface. 
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Figure 3.18 Comparison of Steady and Unsteady Velocity Distributions 
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Smoothness should be ensured in the connecting areas. A compressibility correction may 


also make a difference. 


D. USER’S GUIDE TO BL2D 


1. Output from the Program PANEL 

Three files are output from the panel code described in the previous chapter. 
The file vel.dat is simply a printout of x/c and velocity relative freestream data, mainly 
for plotting purposes if desired. Even though, during computation in the program, panels 
on the lower surface downstream of the stagnation point were taken to have negative 
"directional" velocities, the output file correctly displays the positive "physical" 
velocities. The file cp.dat lists x/c and pressure distribution data. The other output file, 
called bl2d.dat, is generated to be compatible with the program bl2d as input. An 


example input file for the modified program is included in Appendix A. 


2. Input Description 

The first line of bl2d.dat consists of Reynolds number, transition location 
(x/c) on the upper surface, and transition location on the lower surface. The next line 
indicates the number of points and the i value of the location of the stagnation point. For 
the modified version, an additional input value is the transition flag. This indicates if the 
run is an initial estimate (0) or the transition values are fairly well determined and valid 
boundary layer calculations are desired. The remaining lines are identical to the velocity 
output file, except that y/c values of the airfoil coordinates are also included. Note that 


these velocities correspond to a particular angle of attack, the one that was specified 
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when panel was run. The actual value of the angle of attack is not listed separately in 
the file. 

The file bl2d.dat may be edited, replacing the originally specified values in 
the first few lines with new values for the desired flow. The Reynolds number is based 
on the appropriate reference length, usually the chord for an airfoil. The transition 
locations may be obtained from experimental data or from a calculation method. In the 
modified version, Michel’s method may be used to initially estimate the location of 
transition onset. If no information is known before using the program, the transition 
values specified for the first run should be large enough to be downstream of the actual 
transition points, yet not so large that the program will not converge. The arbitrary 
values supplied by the program panel are 0.8 for the upper surface and 0.999 for the 
lower surface. If the transition flag at the end of the next line in the input file is set to 
0, a run of the program will show estimates for transition location on the screen but these 
values will not be used for the boundary layer calculations. The input file should be 
edited, replacing the initial downstream transition values with the estimates and changing 
the transition flag to 1. If the program does not converge with these values, it will be 
necessary to experiment with slight deviations in the upper surface value. Since the 
transition value calculated is only an estimate, this experimentation procedure is the rule 
rather than the exception, especially at higher angles of attack. 

For the stagnation location, it is important to not simply take the i value from 
the velocity output which corresponds to the velocity closest to zero. In order to work 


properly, bi2d usually requires the i value for the first point after the stagnation point, 
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where points are numbered in a clockwise direction from the trailing edge. With only 
positive velocities to inspect, it is impossible to tell where this occurs. The program 
panel, however, uses the negative directional velocities in its calculations and 
automatically determines the proper i value to send to the output file bi2d.dat. If 
convergence cannot be obtained by varying the transition location, it may help in some 


cases to decrease the predetermined i value of the stagnation location by one. 


3. Program Operation 
When all necessary values have been changed, rename the input file 

appropriately for reference, such as bi5.dat for an AOA of 5°. Subsequent runs of the 
program panel for other AOA’s will overwrite the output file bl2d.dat. When ready to 
run bl2d, copy the desired input file to bl2d.dat, the required input file name. To run 
the program, type: 

>bld 
Convergence may most easily be observed by plotting output files such as the skin 
friction coefficient. Modifications to various parameters as explained in the previous 
section may facilitate convergence. In some cases, however, such as an unusually shaped 
airfoil, a highly cambered airfoil, or a standard airfoil with a faulty input file, 
convergence may not be attainable. Additionally, all airfoils at a high enough angle of 
attack will cause the program to break down, as the direct boundary layer method cannot 


handle significant separation. 
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Modifications may be made to the resulting plot using gnuplot commands, Appendix B, 


or editing the command file in a separate window. 


5. PC Version 

To increase its accessibility, bI2d was also converted to a version compatible 
with personal computers. Additional programs included on the Boundary Layer 
Analysis disk are a PC version of panel, an airfoil point generation program called 
airfoil, a PC version of gnuplot, and the command plotting file profile. Most of the 
instructions are the same for this version, but there are a few differences. The programs 
panel and bi2d are restricted to 100 or less airfoil points due to array limitations in the 
PC FORTRAN compiler. An automatic rerun feature was incorporated into bl2d for the 
case of initial transition estimates, where the first run internally restarts using the 
calculated estimates for boundary layer calculations. Finally, the velocity profiles are 
output only to a comprehensive file called profilel.dat instead of FORO60.DAT. There 
are no individual velocity profile files. The horizontal line file called FOROS5.DAT in 


the UNIX version is called profile2.dat for the PC version. 
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IV. VISCOUS-INVISCID INTERACTION CODE 


A. THEORY 

The direct boundary layer code, described in the previous chapter, calculates a 
displacement thickness for a prescribed pressure distribution (or equivalently, a velocity 
distribution, for incompressible flow). As the name implies, it is a direct calculation 
involving one pass, thus very little time is required. Another method, known as the 
inverse boundary layer method, calculates a pressure (or velocity) distribution for a 
prescribed displacement thickness. The displacement thickness represents an "effective 
body" as far as the flow is concerned. Iteration is required and the method requires 
considerably more calculation time. The main advantage of the inverse method is the 
ability to calculate through regions of separation. 

The simplest viscous-inviscid interaction method divides the flowfield into an inner 
viscous region where boundary layer calculations are performed and an outer inviscid 
region where potential flow analysis prevails. The solutions are then iterated until they 
match along the dividing line. This method has "weak" interaction because the only 
exchange of information is along the boundary. 

In the strong interaction method, both the pressure and displacement thickness are 
treated as unknowns and are solved simultaneously with successive sweeps over the 


airfoil. The external boundary condition for the boundary layer equation at the outer 


edge of the viscous region is 





U(x) =u, (x) +8u,(x) (4.1) 


where u,°(x) is the inviscid velocity over the airfoil and 6u,(x) is the perturbation due to 


viscous effects, expressed as 


pd *()) 4.2 
bu,(2)=— Fg ltelO8"(0)| — do (4.2) 


where d(u.6*)/do is the blowing velocity. Equations 4.1 and 4.2 comprise the 
interaction, or coupling law. The interaction takes place between x, and x,. The 
integral term is known as the Hilbert integral, in which the displacement effect is 


modeled by source/sink distributions using potential flow theory. 


B. COMPUTER CODE 


1. Overview 
The viscous-inviscid program incorporates a self-contained panel code based 
on the Hess-Smith method, a boundary layer calculation routine, and an interaction 
scheme. The inviscid panel method is similar to the program panel described in Chapter 
2; however, the boundary condition of zero normal velocity on the surface of the body 
is replaced by a blowing velocity determined from the boundary layer calculations. This 
blowing velocity is used to represent the viscous effects on the inviscid flow. The 


boundary layer method is similar to the program bl2d with some modifications that 


account for the presence of the wake and for low Reynolds number flows. The viscous- 





inviscid code used for the investigation was developed by T. Cebeci. Run time on the 
Stardent computer is about five minutes. 
2. Models 
a. Turbulence Model 
As in the direct boundary layer program, the eddy-viscosity formulation 
of Cebeci and Smith is used, with separate formulas for the inner and outer regions. 
This model in this program has the additional features of low Reynolds number effects 


and a wake flow model. The modified equations are expressed as: 
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and the transition length is represented by 


us, 
xy € 
4. 
(3)\Wop, as (4.5) 
C2} v2 +e 
C?=213 flog Re, -4,7323] 

For the wake flow, the eddy-viscosity formulation is 

= (4.6) 
Ep, = (Ey TE ge (€ naen| 205 me 


where (¢,,), specifies the eddy viscosity of the far wake, taken as the maximum of the 


lower and upper wake eddy viscosities: 


Jenin 
(py =0.064 f (u,-u)dy 
* (4.7) 


(J = 0.064 f (u,-u) dy 


Join 
with y,:, the location where u=u,,,. 


b. =‘ Transition Model 
The program uses an empirical formula called Michel’s criterion to 
calculate a first approximation to the transition location on the upper and lower airfoil 
surfaces. It is expressed as a relationship between the Reynolds numbers based on 


momentum thickness and on the x (surface) location of transition: 
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22400 


Re, 


R,=1.174 0+ ~—)Re, (4.8) 


This is the same equation used in the modified version of bd. The program incompbl, 
however, sometimes adjusts the resulting value. If an area of separation with subsequent 
reattachment is calculated, the transition onset is moved to the beginning of the 
separation bubble. The transition location that is printed in the output file always 
corresponds to a nodal point. In many cases, the program will produce a fairly accurate 
result. Sometimes, however, a refinement may be needed. This may be done by 
experimenting with small shifts around the calculated value, as described in Chapter 3 
for the direct boundary layer program. The experimentation may be implemented by 
using the transition specification option, which overrides calculation. Input options are 
described in the User’s Guide section. 

A more advanced technique, suggested by Cebeci, is the e* method, 
which makes use of linear stability theory. A separate stability/transition code 
incorporating this method is run using the output of the viscous/inviscid code (Michel’s 
criterion) as the first estimate. The new value is supplied as input to the first code, this 
time overriding Michel’s criterion. This type of manual iteration continues until 
convergence, usually within three to four cycles. At the present time, determination of 
all required input values and analysis of the output require significant experience and 


judgment, prohibiting a programmed link between the two codes until further refinement 


is accomplished. The e* method was not used in this study. 





As in the program bl2d, the program incompbl uses the Chen-Thyson 
transition range model to calculate the length of the transition region. Equation 4.5 
shows the modified version of this model. The program incorporates an improvement 
for G,,,, the transition length parameter. G,, may be identified by reducing the (3/C’) 
term to the form (1/G,,), with G,,,=C?/3. In the original model, the constant C has a 
recommended value of 60, resulting in G,,=1200. Whereas a value of 1200 may work 
well for large Reynolds numbers, values from 20 to 80 have been shown to be most 
successful in low Reynolds number flows where separation bubbles exist [Ref. 5]. The 
program determines an appropriate value using an empirical correlation formula in the 


form of C’, also shown in Equation 4.5, 
3. Numerical Techniques 


a. Hilbert Integral 
Equation 4.1, containing the Hilbert integral, may be approximated in 


discretized form as 


a 
ux) =ug(x) +O ¢,{u,8*-u,5") 
jet 


where u,“(x) corresponds to the inviscid velocity distribution which contains the 


displacement thickness effect (5*)* and c, is a matrix of interaction coefficients. which are 


functions of geometry only. 





b. FLARE Approximation 
In regions of recirculating flow, such as a separation bubble, numerical 
stability difficulties may be encountered. The FLARE approximation, due to Fliigge- 
Lotz and Reyhner, neglects the longitudinal convective term u(du/dx) in the region of 


negative u velocity. 


C. RESULTS 


Studies were conducted to: 


@ Determine the effect of including the wake in the calculations 


© Investigate the possible occurrence of negative skin friction and to determine its 
significance 


© Compare the results with the direct boundary layer code 


Since the viscous-inviscid method has the abiliry to calculate through regions of 
separation, additional information was anticipated beyond that provided by the direct 


boundary layer method, 


1. Wake Calculations 
Since one. of the input options is for the inclusion of wake calculations, 
investigations were performed to determine its effect. Calculations may be limited to the 
airfoil surface only, or a grid extending into the wake region may also be used. A 
representative sample of the results is shown in Figure 4.1, depicting the skin friction 
coefficient for three angles of attack. At 2°, the main difference is a movement aft of 


the transition point, where the remainder of the curve maintains its original shape. At 
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Figure 4.1 Effect of Wake Calculations 
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6°, the transition point also moves forward; however, the shape of the turbulent section 
is slightly altered. At 10°, the transition point is unchanged. More significant changes 
occur progressing across the top of the airfoil, leading to a delay in the point where the 
skin friction falls below zero. Thus, the primary effect of wake flow is to reduce flow 
separation on the airfoil, especially important at higher angles of attack. This will allow 
calculations to continue when convergence may not have been reached otherwise. These 
results are consistent with those found by Cebeci [Ref. 6]. Therefore, the wake was used 


in all further studies using this code. 
2. Laminar Separation 


a. Comparison with Direct Boundary Layer Method 
The program incompbl was run for the same conditions as the previous 
study with bI2d, a NACA 0012 airfoil with a Reynolds number of 540,000. The angle 
of attack was increased in two degree increments. Figures 4.2 to 4.7 compare the skin 
friction results for both programs. The low angles of attack show excellent agreement 
in the laminar region. As the angie of attack increases, the most notable difference is 
the absence of C; < 0 for the viscous-inviscid method. Since the direct boundary layer 


code failed to converge with separation greater than that produced at 10°, no comparison 


could be performed beyond this point. 





INCOMPBL —— 





Figure 4.3 C, Comparison, NACA 0012, AOA=2°, Re=540,000 
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Figure 4.5 C, Comparison, NACA 0012, AOA=6°, Re=540,000 
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Figure 4.7 C, Comparison, NACA 0012, AOA=10°, Re=540,000 
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b. Investigation of Higher Angles of Attack 
In Figures 4.8 to 4.14, the angle of attack for the viscous-inviscid code 
was increased further in search of separation, indicated by C;<Q. Finer increments were 
used when separation appeared to be imminent. Separation did not occur until 13.3° 
and was gone by 13.4°. A final plot at 14° shows that the remaining trend is an increase 
in the minimum value of C;. The amount of separation at 13.3° was minuscule, as only 


one data point fell below zero with a value of -0.00065. 





Figure 4.8 INCOMPBL: C,, NACA 0012, AOA=12°, Re=540,000 








Figure 4.10 INCOMPBL: C,, NACA 0012, AOA=13.1°, Re=540,000 . 
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- Figure 4.12 INCOMPBL: C,, NACA 0012, AOA=13.3°, Re=540,000 








Figure 4.14 INCOMPBL: C,, NACA 0012, AOA=14°, Re=540,000 





D. USER’S GUIDE TO INCOMPBL 


1. Required Files 
The files required in a directory to run this program are: 
© incompbl*:, a compiled, executable program compatible with the computer being 
used (arbitrary name if starting from source code) 


@ FORO01.DAT (Stardent) or fort.1 (IRIS): first input file containing airfoil related 
data (required name) 


© incompbl.dat: second input file containing flow and calculation information 
(arbitrary name) 


2. UNIX FORTRAN 

There are three ways to get data from external input files into a FORTRAN 
program on a UNIX based system. The first uses an OPEN statement which associates 
a specific input file name with a file number, just as on a personal computer. This 
method is not used in this program. The second method has READ statements using file 
numbers not declared by an OPEN statement. In this case, a default file name such as 
FORO01.DAT is assumed for a READ statement referencing a file number of one, for 
example. The program uses this method for the first input file (hence the required name) 
as well as many output files, which the user can modify or add to in the source code 
before compilation. The third method, which may be used only once in a program, is 


to specify an input file on the command line at execution time. The file can have any 


name chosen by the user as long as that file has data in the correct format expected by 





the program. READ statements using a file number of five assume interactive input 


from the keyboard during execution or an input file name in the command line. 


3. Starting from the Source Code 

If the compiled version is not available or it is not certain which source code 
file corresponds to the compiled file that exists, the source code is the best starting point. 
Otherwise, proceed with execution procedures in the next section. 

The source code incompbl.f may be obtained from either the Stardent or IRIS 
computer and transferred from one to the other as needed using the ftp utility. Once 
compiled, however, the program will be computer specific, even if the executable 
versions have the same name. 

Modifications may be made to the source code first if desired. For example, 
a WRITE (32,*) statement could be inserted to output specific data to a file with a 
default name of FOR032.DAT (Stardent) or fort.32 (IRIS). This method is often used 
to quickly obtain files of x/c vs. velocity, pressure, or other parameters for plotting. 

To compile the FORTRAN source code for the Stardent, type: 

> fe -O2 -o incompb! incompbl.f 
The term fe is for FORTRAN Compile. The term -O2 permits vector optimization. The 
program will run without it, but a warning notice will be issued after compilation. Note 
the use of the letter "O", not the number zero, "0". DO NOT use the option -O3 
instead. This option is for parallel processing, a feature not currently incorporated on 
the NPS Stardent. The program may appear to successfully compile and run, but there 


will usually be errors in the output. The term following the -o is the name of the output 
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executable program. Any name may be used. The program name will be displayed with 
an asterisk (*) following it in the directory listing to indicate that it is an executable 
program. Finally, the source code or codes are listed in order, only one being used in 
this case. To compile on the IRIS, type 
> {77 -O3 -o incompbl incompbl.f 

The optimization levels have different meanings on the IRIS and -O3 is the correct 
parameter. 

During the compilation process, a file incompbl.o, called an object file, is 
produced appears in the directory listing. This file is not needed in this application and 


may be deleted. 


4. Input File Editing 

The first input file pertains to the panel method part of the program and is 
called FORO01.DAT (Stardent) or fort.1 (IRIS). A sample file is included in Appendix 
A. The first line is simply a number telling how many of the following lines are for 
comments, The next few lines contain the comments, such as the type airfoil being 
analyzed, the date of the test, or any other information useful to the user. The next 
group of data consists of the angle of attack (ALPI), the x/c location of the pivot about 
which the airfoil rotates to a new angle of attack (PIVOT), and the number of panels 
defining the lower and upper surfaces (NLOWER and NUPPER, respectively). Finally, 
the x/c and y/c coordinates are listed in separate blocks, with the order starting at or near 
the trailing edge, proceeding across the lower surface, then the upper surface, and ending 


at or near the trailing edge. The number of points will be one higher than the number 
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of panels, even if the first and last points coincide. The trailing edge point simply is 
listed twice in this case. The only part of this file that is likely to get changed on a 
routine basis is the angle of attack, as the other values are usually fixed for a given 
airfoil. 

The other input file pertains to the boundary layer part of the code. This file 
is normally called incompbl.dat, although the user may give it any name and use this 
name on the command line at the time of execution. A sample input file is included in 
Appendix A. IWAKE is the viscous wake flow flag. A zero indicates that these effects 
will not be included, while a one indicates that they will be included. NXT dictates the 
number of chordwise stations on the body. NW sets the number of chordwise stations 
in the wake. ITREND refers to the number of calculation cycles, where 20 is a good 
starting number, and 30 or 40 may be needed. ITR(1) is a flag for the transition location 
specification method for the upper surface. A zero will activate a calculation using 
Michel’s method, and a one indicates that the location will be provided as part of the 
input. ITR(2), for the lower surface, should be zero. ISWPMX< is the number of sweeps 
in each cycle. A cycle corresponds to the calculation of inviscid and viscous flow 
equations. One sweep is usually sufficient but, in some cases, it may be necessary to use 
2 or 3 sweeps in one cycle. RL is the Reynolds number based on chord length. 
XCTR(1) is the x/c value for the transition location on the upper surface. This value is 
only used if ITR(1)=1; otherwise, it will be ignored. IP is the print flag, which should 


normally be set to one to obtain output. This screen output can be redirected to a file 
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for later use by including the proper option on the command line at the time of 


execution. If IP=0, the standard output will not be generated. 


5. Program Execution 
After ensuring the appropriate files are properly edited and are present in the 

directory, type: 

>incompbl <incompbl.dat > incompbl.out 
Several new files will appear in the directory after running the program. Any write 
statements incorporated into the program using the default numbering format will produce 
files such as FORO32.DAT (Stardent) for fort.32 (IRIS). The output file incompbl.out 
will have a comprehensive summary of the input data, and output data such as C,, Cp, 
C_, Cy, and boundary layer properties, including skin friction coefficient and 


displacement thickness. 
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Vv. EXPERIMENT 


A. BACKGROUND 

The experimental work described in this chapter was performed as part of a 
continuing series of investigations by M. S. Chandrasekhara and L. W. Carr in the 
Compressible Dynamic Stall Facility (CDSF) of the Fluid Mechanics Laboratory (FML) 
at NASA Ames Research Center. The steady flow density field around a NACA 0012 
airfoil at a Reynolds number of 540,000 was photographed using the technique of Point 
Diffraction Interferometry [Ref. 7]. 

A sample photograph, called an interferogram, is shown in Figure 5.1. The bright 
and dark areas emanating from the airfoil are called fringes, and each one represents a 
line of constant density. The stagnation point may be easily identified as the center of 
the smallest fringe on the lower surface (for a positive angle of attack). The flow 
accelerates around the leading edge of the airfoil. The example shown also reveals the 
presence of a laminar separation bubble just aft of the leading edge, distinguished by a 
characteristic fringe pattern. 

For a given freestream Mach number, the fringes may also be correlated with 
particular Mach numbers and pressures in isentropic flow. For a standard interferometer 


and two dimensional flow, the path length difference APL due to density changes can be 


related to the fringe number e: 








Reyes 


ure 5.1 Interferogram of NACA 0012 Airfoil, AOA= 0° 





Fig 
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APL L 
€=—— = (4-4 — 5.1 
acs (6.1) 


where n is the refractive index of the signal beam, n,.- refers to the reference beam, Ag 
is the wave length of light used, and L is the test section span. An integer value for e 
results in a bright fringe, while an integer plus one half corresponds to a dark fringe. 
Using the Gladstone-Dale equation [Ref. 8] and the perfect gas equation, Equation 5.1 


teduces to 





Ao . €Ro 
= ps 9). 5.2 
P-P., mri rae (5.2) 


where p is the density of the fringe, p, refers to freestream conditions, and A is a 
constant determined from the experimental parameters. With \)»=532 nm, L=25 cm, 
(no-1)=2.733 x 10% , and total or stagnation density p>=1.21 kg/m’, A=0.009421 


kg/m’®. Dividing by po, 





ib Ph Be (5.3) 
Po Po Po 
Using the relationship 
ae 
fats sae ul (6.4) 
Po 2 


with p=p,, and M=M,,, p../po iS a function of freestream Mach number only for y=1.4. 


For the present case of M=0.3, the term p,./pp is constant at 0.956. Equation 5.3 may 


now be written as 





P. =0,956 +0.007786€ (5.5) 
Po 


Thus, quantitative density measurements are available knowing only the fringe number. 
The fringe numbers are determined by identifying the fringe number of the stagnation 
fringe as the highest fringe number, and counting down around the leading edge and 
upper surface of the airfoil. With p/p)=1 in Equation 5.5, «+6. Knowing the densities, 
pressures may be calculated in a straightforward manner: 


PD. al (5.6) 
Po \ Po 


Mach numbers for each fringe may be calculated using Equation 5.4. 
The photographs for various angles of attack are identified using particular settings 
on the equipment. A correlation between settings and angles of attack is presented in 


Table 5.1. 


Table 5.1 CORRELATION OF EXPERIMENTAL SETTINGS AND AOA 


DEGHEES” COUNT | DEGHEES COUNT [DEGREES COUNT] 








108 14.5 165 
15 171 
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B. IMAGE PROCESSING 


1. Scanning 
The processing task began with a set of experimental interferograms at angles 

of attack ranging from 0° to 10°. Each was scanned into an IRIS computer using the 
program pixscan at the Numerical Aerodynamic Simulation Facility at NASA Ames 
Research Center. Options for contrast enhancement (gamma correction) and grayscale 
were used. A photograph of the airfoil in no-flow conditions was also scanned. All of 
the photographs were carefully placed on the scanner in a position such that the top 
surfaces of the two bottom reference triangles made a horizontal line even with the edge 
of the scanning window. The resulting scanned images were transferred via ftp to the 
IRIS computer at the Naval Postgraduate School for further processing. 

2. Editing 

Using a program called pixedit, the original flowfield images were first overlaid 
with the airfoil image. This permitted a more defined surface for analysis when the 
diffraction caused by the interferometry technique distorted the airfoil surface. The 
images were then cropped closer to the airfoil. This procedure reduced the number of 
pixels in the image, which was necessary for use with the digitizing program. A small 


section near the leading edge was also selected for enlargement, thus giving better detail 


of the laminar separation bubble region. 





In anticipation of the digitization process, reference marks at known coordinates 
were added to the two images. The original photographs provided the basis for 
coordinate system definition. The three dark triangles are located as shown in Figure 
5.2. The distance between horizontal or vertical faces of the triangles is 0.2 of the 
chord. The vertical faces of the two triangles on 
the right are located at 0.25 chord. Thus, the 


coordinates of the point at the right angle of the 


PARALLEL TO. Maca FACES 
biworte 28% CHORD LOC: 


bottom right triangle are (0.25, -0.1), A 


AIPEONCENTERUNE 


reference mark is needed in the top left corner une raraes ronorzoNraL races 
s RRA 30 CHORD Lint 


for use with the digitization program, but none is 
\ — —_—- | / 
CHORD 


provided. Therefore, a mark was constructed 





with several applications of the snap new image Figure 5.2 Reference Triangles 

_ feature of pixedit. Very thin horizontal and vertical lines were initially saved as separate 
images. For each interferogram, these images were read in and placed to be even with 
existing reference marks of the airfoil leading edge (vertical line) and the bottom edge 
of the upper triangle (horizontal line). The small area in the top left corner where the 
lines crossed was saved as another image and the long lines were deleted. The 
coordinates of the cross mark are (0.0, 0.1). Finally, the complete image was saved for 
further processing. A similar procedure was used to make a cross mark at (0.05, 0.0) 


for the lower right reference on the enlarged leading edge images. 





3. Fringe Tracing 
Both of the new images, the density flowfield and the enlarged leading edge, 
were digitized using the program DigiCurv. The appropriate lower right and upper left 
coordinates were provided and the corresponding points on the image were selected with 
the mouse. With the coordinate system thus initialized, each curve was digitized 
separately. The centerline of the dark fringe was used for digitization. A representative 
number of points were chosen for each fringe using the left button of the mouse, always 
starting with the point where the fringe intersected the airfoil. Depressing the right 
button activated a menu with a curve fit option. The program calculated a best fit curve 
to the chosen points using a spline routine and output up to 30 (default value) new points 
to describe the curve. The fitted curve was also displayed on the screen for acceptance. 
In most cases, the computed curve had outstanding agreement with the centerline of the 
dark fringe, even when a sharp corner was involved near the end of a separation bubble. 
If the displayed curve needed refinement, the backup option removed points one by one, 
more closely spaced points were selected, and the curve was refit. After all dark fringes 
were digitized for an image, the data point output file was saved for plotting. 
The digitized data can be viewed using any plotting routine compatible with 
sequential data separated by blank lines. The program xyplot does not recognize blank 
lines and is not a good choice. The program gnuplot will properly accept data in this 


format. The digitized plot corresponding to the interferogram in Figure 5.1 is shown in 


Figure 5.3. 
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Figure 5.3 Digitized Interferogram of NACA 0012 Airfoil, AOA=10° 
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4. Printing 
To obtain graphical output on a postscript laser printer, the files must be in 

postscript format. The digitized plots are in this format and can be easily printed using 
the following commands in the Advanced Computation Laboratory at the Naval 
Postgraduate School. From the Stardent, type 

>rglp filename 
From the IRIS, type 

> Ip -dlaser filename 
The filenames may have the suffix ps to indicate that they are postscript files, but this 
is for reference only and is not required by the printer. 

The images output from pixedit are not directly compatible with the postscript 
printer. Using the IRIS computer, they must first be converted to sgi format: 

> pix2sgi infilename outfilename.sgi 
Again, the suffix is optional. The next step is a conversion to postscript format: 

> tops infilename.sgi -p 98.0 > outfilename.ps 
All suffixes are optional. The -p 98.0 option preserves the proportions of the original 
image; otherwise, it may be distorted when stretched to fill up the printed page. 
The > symbol redirects the screen output to the specified output file. There are a 
number of other options available, including size and color. For a full description of 
options, type 


> man tops 





to obtain the on-line help manual entry on this conversion program. 
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Cc. ANALYSIS 


The primary information desired from the interferograms was the location of the 
start of the laminar separation bubble, when it existed. Figure 5.4 shows an enlarged 
image at 6°, where a bubble is first formed. Using the digitized data file of this image, 
precise x/c locations of the intersections of the fringes with the airfoil surface may be 
determined, as is the first point digitized for each fringe. Starting with the stagnation 
fringe, the first dark fringe always has a number of 5.5 (for the present case of 
M,,=0.3). Subsequent dark fringes have numbers 4.5, 3.5, etc. Negative fringe 
numbers are possible. Converting pressures to pressure coefficients (C,), a table of 
fringe numbers and corresponding C,’s may be produced and used for all interferograms 
with the same freestream Mach number. The program pres (Appendix A) was written 
to process the data as described and output a plotting file of x/c vs. C,. Figure 5.5 
shows the digitized image at 6° and the corresponding C, plot. The most important 
feature is the pressure plateau, which reveals the existence and location of a separation 
bubble. The bubble starts at a chord location of 0.015. Interferograms for other angles 
of attack were analyzed in the same manner. Figure 5.6 shows a comprehensive plot for 


an angle of attack range from 6° to 10°. As the angle of attack increases, the length of 


the bubble increases, and the starting location moves slightly forward. 








Figure 5.4 Enlarged Leading Edge Int erogram, AO ‘ 
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Figure 5.5 Digitized Interferogram and C, Plot, AOA=6° 
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Figure 5.6 C, Plots Showing Separation Bubble, AOA 
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D. COMPARISON OF RESULTS TO COMPUTATION 

Table 5.2 compares the experimental results to those obtained by bl2d and 
incompbl. The direct boundary layer code shows a very small separation bubble at 4°; 
however, none is present in the experiment. At 6°, the bubble’s existence is correctly 
detected, but the computational location is slightly aft of the experimental value. The 
higher angles of attack show excellent agreement, with differences of only 0.3% of 
chord. In contrast, the viscous-inviscid code failed to predict any separation until 13.3°. 
Since this is past the steady stall angle of 12.4°, there is no experimental bubble to 
compare with. Even if stall had not occurred, the trend clearly indicates that the location 


would be in great error as well. 


Table 5.2 COMPARISON OF BUBBLE START LOCATIONS 
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VI. NAVIER-STOKES CODE 


A. OVERVIEW 

A time-averaged Navier-Stokes (NS) code called ns2 was used for a final study. 
This method has the advantage of including compressibility effects. While at a Mach 
number of 0.3 the effects are small, it is just on the border of the region that is usually 
considered acceptable for the assumption of incompressibility. A disadvantage of this 
method, as mentioned in the transition discussion of Chapter 2, is the lack of a transition 
model. Turbulent flow is assumed throughout the flowfield. The code also takes three 
to four hours to run on a Stardent computer. Nevertheless, Navier-Stokes codes are 
often regarded as the best computation method currently available and warrant 
consideration. The details of the time-averaged NS equations, their derivation, and their 
discretization are well documented elsewhere [Ref. 9] and will not be reviewed here. 
The particular code used for this investigation was developed by J. A. Ekaterinaris of the 


Navy-NASA Joint Institute of Aeronautics. 


B. RESULTS 
The code was run for a NACA 0012 airfoil, with a Mach number of 0.3 and a 
Reynolds number of 540,000. A 161x64 C-type grid with a very fine distribution normal 


to the surface in the viscous region was used. The grid is shown in Figure 6.1. The 


Baldwin-Lomax model was used for turbulence modeling. The program was run to 2000 
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Figure 6.1 161x64 Viscous Grid for NACA 0012 Airfoil 
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iterations, where the residuals had dropped two orders of magnitude. Figures 6.2 to 6.4 
show comparisons of the computational density fields to the experimental interferograms. 
At 0°, there is no separation bubble and the agreement appears to be good. At 6°, the 
experimental bubble first appears at x/-=0.015. All of the computational density curves 
emanate from the leading edge. This is consistent with the lack of a transition model. 
At 10°, both methods show a large bubble. The NS bubble is much further aft. There 
is also a small extra bubble near the leading edge. It appears that the "state of the art" 
method is not always best for representing reality, particularly for low Reynolds number 
flows. The author is, however, currently working on a version of the program which 
will incorporate transition calculations and preliminary work indicates that results will 


be much better. 
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Figure 6.2 Comparison of NS and Experiment, AOA=0° 
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Figure 6.3 Comparison of NS and Experiment, AOA=6° . 
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VII. CONCLUSIONS 


In the computational investigations of a NACA 0012 airfoil in low Reynolds 
number flows, several important discoveries have been made. First, even though Navier- 
Stokes codes are the most advanced computational method currently available, they are 
not always appropriate. Specifically, a transition model is necessary for the case of low 
Reynolds number flows. Even if accurate results may be obtained, the method is not 
practical for many applications due to its high cost in time and money. The viscous- 
inviscid method seems to offer very advanced calculations at a very inexpensive price. 
The problem with this code used alone, however, is that it simply does not give correct 
results for separation bubbles in low Reynolds number flows. A compressibility 
correction in a future version may alleviate the problem. Using a stability/transition 
method in conjunction with the code will certainly provide more refined transition 
estimates, which may influence separation bubble results. Finally, the relatively simple 
and often overlooked direct boundary layer method can provide meaningful information 
about separation bubbles in low Reynolds number flows. Given that the code is also very 
fast, it may be used efficiently in the design stages and quality assurance checking of 
many aeronautical applications. 

Successful or not, all current computational methods have several important 


limitations. All are dependent on empirical models for transition onset, transition length, 


and turbulence. The models are often formulated for specific parameter ranges outside 








of which agreement is poor. Until there is sufficient computing power available to solve 
. the full Navier-Stokes codes, it is imperative to check the applicability of the models and 


the assumptions of the equations on which a method is based before using a program and 


counting on the results to be reliable. 
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APPENDIX A: COMPUTER PROGRAMS AND FILES 


The following programs and input/output files are presented in logical order of use. The 
input and output files are examples only. The input values may be modified as described 
in the appropriate User’s Guide section. The program incompbli comprises 90 pages in 
length and was not modified other than inserting WRITE statements to extract data for 


plotting; thus, the source code is not included. 
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1.000000 
€.999013 
0.996057 
0.993144 
€.984292 
0.975528 
0.964888 
0.952414 
0.938153 
0.922164 
0.904508 
0.885257 
0.864484 
0.842274 
0.818712 
0.793893 
0.767913 
0.740877 
0.712890 
0.684062 
0.654508 
9.624345 
0.593691 
C.562667 
0.831395 
0.500000 
0.468605 
©. 437333 
0.406309 
0.375655 
0.34549) 
0.315938 
9.287110 
0.259123 
0.23208? 
0.206107 
C.181288 
0.157726 
0.135516 
6.114743 
C.095491 
0.077836 
0.061847 
C.047586 
0.035112 
0.024472 
0.015708 
0.c088aseé 
0.003943 
0.000987 
0.000000 
0.000987 
0.003943 
0.008856 
0.015708 
0.024472 
0.035112 
0.047586 
0.061847 
0.077836 
€.995493 
0.114743 
0.135816 
0.157726 
0.181288 
0.206107 
0.232087 
0.259123 
0.28711¢ 
0.315938 
0.345492 





0.000000 
-0.000143 
~0.000562 
-0.001258 
-0.002222 
-0.003443 
-0.004909 
0.006603 
-0.008510 
-0.010610 
-0.012883 
~0.015310 
-0.017868 
-0.020535 
-0.023291 
~0,026111 
-0.028974 
~0.031856 


~0.040378 
-0.043094 
-0.045705 
-048182 
-9.050499 
-0.052625 
-0.054534 
-0.056195 
-0.057581 
-G.058666 
-0.059424 
-0.059834 
-0.059876 
0.059535 
-0.058799 
-0.0576613 
-0,056119 
-0.054176 
-0.051839 
~0.049121 
-0.046037 
0.042608 
-0.038854 
-0.034800 
~0.030471 
70.025893 
-€.021088 
-0.016072 
-€.010884 
-0.005521 
2.000000 
0.005521 
¢.010884 
0.016078 
0.021082 
2.025893 
0.030471 
0.034800 
0.038854 
0.042608 
0.046037 
0.049121 
0.051839 
0.054376 
0.056119 
0.057661 
0.058799 
0.059535 
0.059876 
0.059834 
0.059424 





Input file: points.dat 


0.375655 
0.406309 
€.437333 
0.468605 
C.S00000 
0.531395 
0.562667 
0,593691 
0.624345 
0.654509 
0.684062 
0.71289¢ 
0.740877 
0.767913 
0.793893 
0.818712 
0.842274 
0.864484 
0.885257 
0.904509 
0.922164 
0.938153 
0.952414 
0.964888 
0.975528 
0.984292 
0.991144 
0.996057 
0.999023 
1.000000 


0.058666 
0.05758) 
0.056195 
0.954534 
0.052625 
0.050499 
0.048182 
0.04570S 
0.043094 
0.040378 
0.037582 
0.034733 
0.031856 
0.028974 
6.026112 
0.023291 
0.020535 
0.017868 
0.015310 
0.022883 
0.010610 
0.008510 
0.006603 
0.004909 
0,003443 
0.002222 
9.001258 
0.000562 
0.000141 
0.000000 
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Source code: panel.f 


PROGRAM FANEL 


AUTHOR: L. M. NOWAK 





MFTHOD: 


* NUMBER OF NODES ON AIRFOIL SURFACE: 


PI=ACOS(- 


6 NOV 93 


LIM: Arrays currently dimension 
Input data file points.dat 


modified: MAY, AUG 1992 
CALCULATE THE VELOCITIES ON AN AIRFOIL USING A PANEL METHOD. 


ed for maximum of N=200 panels 


will have N+) points 


Output velocities are referenced to freestream, ie. V/vinf 


PANEL. VORTICITY STRENGTH 
BOUNDARY CONDITIONS INCLUDE 


FLOWFIELD CONSISTS OF THREE SIMPLER FLOWS: 
AND VORTICITY. SOURCE DISTRIBUTIONS q{j) VARY FROM PANEL To 
GAMMA IS THE SAME FOR ALL PANELS. 
FLOW TANGENCY AT CONTROL POINTS AND 


KUTTA CONDITION FOR FIRST AND LAST PANELS, - INFLUENCE 


COEFFICIENTS COMBINED TO Fo. 
OF n+1 EQUATIONS, n+1 UNKNOWNS (q(1)...q(n), GAMMA). 


AT CONTROL POINTS EVALUATED FROM q{j) AND GAMMA, 


REAL % (12202), ¥(1:202), XM) 


3202}, YM(1:202), 


At (1:202,1:202), Bt (1:202,1:202), 

(25202, 1:202),b(1:202), 

922202), ve 115202) ,ALPRA, Vv, VEC (11202), 

PI,GAMMA, THETA (1:202) , NUM, DEN, 

R(1:202, 12202}, BETA(1:202,1:202) ,NUM1, DENI, NUM2, DEN2, 


AAUG (2 2202,1:202),An(1 


2202,1:202},Bn(1:202,1:202) 


PRINT*,‘ INPUT NC, OF PANELS {1 less than flines in points.dat):’ 


READ *,N 










CONST FUL 


print *,’ INPUT REYNOLDS NUMSER:¢ 


print *, "ENTER C 
PRINT -* 5+ : 
READ *, IANS 


READ *,RE 






IF(IANS,EC.1) THEN 


PRINT *,*INPUT X/C TRANSITION 
READ *, TRANSUPPER 
PRINT *|¢INSUT X/C TRANSITION 
READ *, TRANSLOWER 


ELSE 
***These are arbitrary valves intended 
actual transition points, for use 


ENDIF 
WRITE (90,50) 


TRANSUPPER= .& 
TRANSLOWER=.999 


.dat’, STATY, 
ot.dat’, STAT 


RL, TRANSUPPER, TRANSLOWER 












* UNKNCHN’ } 
UNKNOWN® } 





fi, lat’, STATUS=" UNKNOWN’ ) 
OPEN [UNIT=90, FILEs*bi2d.dat’, STATUS=" UNKNOWN! : 






FF TRANSITION LOCATIONS UNKNOWN’ 
IF TRANSITION LOCATIONS KNOWN: ' 





LOCATION FOR UPPER SURFACE:* 


LOCATION FOR LOWER SURFACE:’ 


to be downstream cf the 
with Michel’s criterion in aL2D 







WRITE (94,50) RL, TRANSUPPER, TRANSLOWER 
FORMAT (F10.0,F10.4,F10.4} 


PRINT *,* INPUT ANGLE OF ATTACK IN DEGREES:’ 


READ *,ALPRA 
ALPHA=ALPHA*PI/180.0 


DO 30 1=1,Ns1 


CONTINUE 


READ (88,25) X{1},¥(1) 
FORMAT {(2(F8.6,2x))} 












FREESTREAM, SOURCE, 


RM NEW COEFFICIENTS IN LINEAR SYSTEM 





VELOCITIES 












































*This section defines the influence coefficients: 


DO I=1,N 
XM(1) 20.5" (X (1) +X (14+2)) 
YM(T)=0,5*(y (1) +y (Tol) } 
RUT, 1.5 (XM (2) XK C1P) 82.4 YM CTP -¥ (4) ) 82, 
DO J=1,N 
NUM=Y (J+) ~¥ (J) 
DEN=X (J+1} -X (J) 
THETA (J) =ATAN2 (NUM, DEN} 
NUM1=YM(I)-¥ (J+1) 
DEN1=XM(I)-X(J+1} 
NUM2=YM(1)~¥ (J) 
DEN2=XM({I)-X{J} 
BETA(T, J) =ATAN2 { (NUM) *DEN2-DEN1*NUM2) , (DEN1*DEN2+NUM! *NUM2) 
R(I,d+h XM CLD -X (D419) 882.4 CYM (1) -¥ (S41) AZ, 
THETADIF=THETA (I) ~THETA (J) 
IF (1.EQ.3} 
THEN 
An(I,J)=0.5 
Bn{t,J)=0.0 





ELSE 
An(I,d)=(1/ (2*21) ) * (SIN (TRETADIF) *ALOG(R(Z,J#1) /R(T,I)} 
* .$+COS(THETADIF) *BETA(I,J)) 
Bn{I, J}=(1/(2*PI)) * (COS (THETADIF) *ALOG(R(I,J+1)/R(1,5}) 
* S~SIN(THETADIF) *BETA(I,J)} 
ENO IF 
At (1,J}=-Bn(I,d} 
Bt (I, J}=An(1,J) 
END DO 
END DO 


* Matrix coefficients of linear system defined (a’s and b's): 


a{N+1,N+1)=0.0 
DO I=1,N 
a(I,N+1)=0.0 
DO J=1,N 
a(l,J)=An(i, J) 
a(Z,Nel}ea(I,Ne1) +an(1,7 
END DO 
b{I) =-1.0*°SIN(ALPHA-THETA(I)) 
a (N41, I)=Ac 1,1) +At (N, 1) 
a(N41,Ns])=a(Ns1,N+1}+Bt (2,7) 4BL (N, I) 
END BO 
bIN+1} =-1.04 (COS {ALPHA~THETA (1) ) +COS (ALPRA-THETA(N} } 


* Define augmented matrix for {input to linear solver subroutine GAUSS 


DO I=1,Ns1 
DO Je1,Ne1 
ARUGII, 3) =a (I,J) 
END DO 
AAUG{I, N+2) =b(T) 
END DO 


CALL GAUSS (N+1,AAUG 
* Define source and vorticity strengths: 
DO I=1,N 
q(I} =AAUG(T,N+2} 
END DO 


GAMMA=AAUG (N+1,N+2) 


* Calculate velocity on each panel at control point 


NSTAGFLAG=0 
ISTAG=0 
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Page 3 


END DO 
Vt (1) vt (I) +COS (ALPHA-THETA (1) ) 
Cp21.d-vt (1) **2 
WRITE (25,*) XM(I),Ve(1) 
IF ((Vt(I}.GT.0) .AND. (NSTAGFLAG.EQ.0)) THEN 
ISTAG=I 
NSTAGFLAG=1 
ENDIF 
Ie (Vt(1).LT.0} Ve(T)s-vett) 
WRITE (89,45) XM(I), Ve (I) 
WRITE (91,45) XMII),-Cp 
END DO 


FORMAT (2(F10.5,2x)) 
FORMAT (3(F10.5)) 
FORMAT (315) 


WRITE (90,49) N,ISTAG, IANS 
DO I=1,N 

WRITE (90,48) XM{I),YM(I),Vt(L) 
END DC 


print *,’CALCULATIONS COMPLETE’ 
PRINT *,/OUTPUT FILES ARE veloc.dat, cp.dat, bi?d.dat' 


END 


PerewererreverTTTrerrrrererri rt tt iri ii tits tee 


* Gauss elimination procedure obtained from Numerica! Methods text ch.6 


SUBROUTINE GAUSSIN, 2) 
INTEGER PV 
REAL 2(1:202,1:203),£ 


£=1.0 
TF (1,0+6.GT.1.0) THEN 
£=€/2,.0 
GOTO 10 
END iF 
EsE*2 
EPS2=2*E 
FRINT *,7 MACHINE EPSTLON=’,£ 


DET= 1 


DO 101C 3=1,N-1 
PV=I 


DO 1020 J=I+1,N 
IF (ABS(Z(PV,I))  . LT. ABS(2(J,1))) PV=T 
CONTINUE 


IF (PV.EQ.1) GOTO 1050 


DO 1040 JC=1,N+1 
TM=Z(1, JC} 
2(1,J5C) =2 (PV, JC) 
Z(PV, JC} =T™M 

CONTINUE 


OET=-1*DET 
IF (2(1,1}.€@.0) THEN 


GOTO 1200 
END IF 
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aaa 


1075 


1060 


1120 
1130 


1200 


DO 1060 JR=I+1, N 
IF (2(3R, 1) .NE.O) THEN 
R=Z(0R, 1/211, 0) 
pO 1075 KC=i+1,Nel 
TEME=2 (IR, KC) 


249R, KC) =Z {IR, KC} -R*Z{1, KC) 
IF (ASS (2 (IR, KC)) w LT EPS2*T: 


EMP} 2(IR, KC}=0.0 


t-- if the result of subtraction is smaller than 








CONTINUE 
END IF 
CONTINUE 
CONTINUE 


Do 1084 I=1,N 
DET=DET*Z (1,1) 
CONTINUE 


PRINT * 
PRINT *,’DETERMINANT = *,DET 
PRINT * 


IF (2(N,N).£0.0) GOTO 1200 
ZIN, Nel) =Z(N, N21) (ZN, 


po 1130 NV=N-1,1,71 
VAs Z (NV, N+3) 
DO 1120 K=Nv+t,N 
VA=VA-2 (NV, KD TZ (K, Ne 2) 
CONTINUE 
2 (NV, N41) =VA/Z (NV, NV) 
CONTINUE 
RETURN 


PRINT *,’MATRIX IS SINGULAR’ 
PRINT sveqee, 1, /20, 08 2D 
STOP 

END 
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- 2 times machine epsilon times the original 
ta- value, it is set to Zero. 












0.19365 
0.16945 
14660 
12510 
1 
0 


co 
RABO 
Soa 
ABRO 
Boow 


= 
eS 
w 
3° 


ooaeooeooo0e°o 
oe0°o 

Senn 

ANCS 

BNO 

retv ree 


oo 
20 
oo 
ON 
ang 
wo 


0.00045 
0.00240 
0.00635 
0.02225 
2005 
.02975 
0.04130 
0.05465 
0.06989 
0.08660 
0.10505 
0.12510 
0.14660 
0.16945 
0.19365 
0.21905 
0.24555 
0.27310 
0.30150 








0.05510 
0.05685 
9.05815 
0.05910 
0.05965 
0.05980 


0.762 


0.75382 
0,90956 
0.84584 
0.87525 
0.89440 
0.91199 
0.92807 
0.94058 
0.95235 
0.96041 
0.97159 
0.98018 
0.98850 
0.99553 
1.00410 
1.01099 
1.01700 
1.02456 
4.03023 
1.03612 
1.04173 
1.04948 
1.05502 
1.05902 
1.06629 
1.07290 
1.07664 
1.08126 
1.08640 
1.09167 
1.09584 
1.09756 
1.10065 
1.10082 


1121308 


Output/input file: bl2d.dat 


0,33C65 
0.36050 
0.39095 
0.42180 
0.45295 
0.48430 
0.51565 
0.54695 
0.57810 
0.60895 
0.63940 
0.66925 
0.69840 
0.72680 
0.75435 
0.78085 
0.80625 
0.83045 
0.85330 
0.87480 





0.99950 
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0,059€0 
0.05900 
0.05805 
0.05680 
0.05530 
0.05355 
0.05250 
0.04925 
0.04690 
0.04435 
0.04365 
0.03890 
0.03610 
0.03325 
0.03035 
0.02750 
0,02465 
0.02185 
0.01915 
0.01655 
0.01405 
0.01170 
0.00955 
0.00755 
0.00575 
0.00415 
0.00280 
0.00170 
0,00085 
0.00030 
0.00005 


0.75382 








Source code: bi2d.f 
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POeVESUSSUSCSVESTCSS OSS OS ITT C tite ttt rstserrs ereerrsy 


Modifications: L. M. NOWAK 


ver. 2 
16 July 1992: 


* * 
* . 
. Added write statements (30-40) in the “DO 175" loop to * 
. output boundary layer profiles for plotting . 
§ {100 panel airfoil only) * 
. 20 July 1992: Added write statements (20-21) to ovtput * 
e CF(skin friction} and DLS{delta star) for plotting ® 
be ver. 3 * 
* 31 Aug 1992: Added calculation for onset of transition * 
* based on Michel’s criterion, added input ITRANS * 
* ver. 4 ‘ 
* 3 Sept 1992: Modified boundary layer profile output = 
* to be compatible with airfoil of any number. points * 
+ ver. S = 
* 3 Sept 1992: Redimensioned al! arrays to accept * 
* airfoil up to 200 panels . 
RERARAAERAM ROSE AERDER ARSE AR ASA DAE SD AEH SERERERERERERANEREERRE REED 

SUBROUTINE BL 

COMMON /BLC2/ NX,NXT,NP,NPT,NTR, IT, ISF 

COMMON X{200} , VE (200) ,P1 (200) , 2 (200) , GMTR(200) 

COMMON £TA(203),DETA(201) ,A(201) 

COMMON F(201,2),U(201,2)},V(201,2) ,B(201,2) 

COMMON DELF (201), DELU(203) , DELV(201) 


NP = ALOG((ETAE/DETA(1)) * (VGF-1.0) +1,0) /ALOG(VGP)+1.001 
ETA(1) = 0.0 
DO 10 J=2,NPT 
ETA(J) = ETA(J-1) + DETA(J-1) 
DETA(J)= VGP*DETA{J-1) 
A(J) = O.S*DETA(2-1) 
CONTINUE 


INITIAL LAMINAR VELOCITY PROFILE 
De 20 J=1,NE 
ETAB = ETA(J)/ETA(NP) 
ETAB2 = ETAB**2 
F(S,2) = 0.25*ETA(NP) *ETAB2*(3.0 - 0.5*ETAB2) 
U(d,2) = 0.5*ETAB* (3.0 - ETAB2) 
V9, 2) 1.5°(1.0 - ETAB2) /ETA(NP) 
B(3,2) 
CONTINUE 


NX 
It 
IGROW 


IT 

IF (IT .GT. ITMAX) GO TO 101 
If (NX .GE. NFR) CALL EDDY 
CALL COEF 

CALL SOLV3 


CHECK FOR CONVERGENCE 

If (NX ,LT. NTR} THEN 
IF(ABS(DELV(1)) .GT. EPSL) GO TO 5 

ELSE 
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=f 
Cc 


30 


Cc 
101 


= 


LE (ABS{DELV(1)/V(1,2)} .GT. EPST) GO TO Ss 


ENDIF 


PROFILES FOR GROWTH 
DO 30 J=NP+1,NPT 
F(J,2) = F(d-1,2) + DETA(J~1) *U{I-1,2) 
U(d,2) utJ-1,2) 
V(J,2) = 0.0 
B(3,2) = B(I-1,2) 
CONTINUE 


CHECK FOR GROWTH 
IF (ABS(V(NP,2)} .GT, 0.0005 .OR. ABS(1.0-U(NP=2,2) /U(NP,2)) 
GT. 0.005) THEN 
NP = NP+2 
IGROW = IGROW+] 
IF (NP .LE. NPT .AND. 
IT = 0 
Go TO S 
ENDIF 
ENDIF 


IGROW .LE. IGROWT) THEN 


CALL OUTPUT 
TF (NX .LT. NXT) GO TO 1 


RETURN 
END 
SUBROUTINE COEF 


COMMON 
COMMON 
COMMON 
COMMON 


/BLC2/ 
/BLC3/ 
/BLCT/ 
/BLCB8/ 


NX, NXT, NP, NPT, NTR, IT, SF 

X (200) , VE (200) ,P1 (200), P2 (200), GMTR{200) 
ETA (201) ,DETA(201) ,A(201) 
F(201,2},U(201, 2) ,V(201, 2) ,B(201, 2) 
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COMMON /BLC9/ 


$1{201) , $2(201), $3(201) ,$4(201} , $5(201) ,56(201), 
+ $7 (201) 58 (201) ,R1{201) ,R2 (201) -R3(201) ,R4(2C? 


) 


P1H = 0.5 * PL(NX) 
IF (NX .£Q. 1) THEN 
CEL = 0.0 
CELH= 0.0 

po 5 J=1,NP 
F(J,1) = 
U(I,1) 
vis,1) 
BtJ,1) 
CONTINUE 
ELSE 
CEL = 0.8 
CELH= 0,5 
ENDIF 


(OX (NX) #X (NX=1)) 7 OXONX) -X(NX= 1) 
CEL 


DO 10C J= 2,NP 
CURRENT STATION 

O.S*(F(S,2) + F(d-1,2)) 
O.5*(U{S,2) + U(J-1,2)) 
O.5* (F(T, 2) *V(S, 2) +P (J-1, 2) *V(I~1,2)) 
O.5*(V(3,2) + V{J-1,2)) 
O.5*{U(T,2)**2 + U{I~-1,2)**2) 
(B(T,2)*V(I,2) = B(J-1,2) *V (3-1, 2)) /DETA(I-1) 


STATION 
O.S*(F(S,1) + F{d~-1,1)) 
C.5*(U(I,1) © U{J-1,1)) 
O.5*(V(I,1} + Vd-1,1)) 
0.5*(U(J,1}**2 © U(J-1,1)**2) 
O.5*(F (3,1) *V(0, 1) 4F (9-1, 1) *V(I-1,1)) 
CDERBV (B(J,1)*V(d,2) = B(J-1, 1) *V(5-1,1)) /OETA(J-1) 
S- COEFFICIENTS 
SL(J) = CELH*(F(J,2) = CFB) + PIR*F(J,2) + 8{J,2) /CETA(J-1) 
S$2(3) = CELK*(F(3-1,2)-CFB) + PLH*F (I-12, 2) -B(3-1,2) /DETA(I-1) 
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$3(J) CELH*{V(J,2) + CYB} + P1H*V¥(J,2) 
34(J) CELH* (V(J-1,2) + CVB) + PlH*V(3-1,2) 
SStJ) ~{CEL+P2 (NX}}*U(J,2) 

S$6(J) ~ (CEL+P2 (NX) ) *U(J-1,2) 


R- COEFFICIENTS 

IF (NX .EQ. 1) THEN 
CRB -P2 {NX} 
R2 (J) CRB - (DERBY * P1{NX)*FVB - P2 (NX) *USB} 

ELSE 
CLE CDERBV * P1(NX=1)*CFYB - P2(NX-1}*CUSB + P2(NX-1) 
CRB -CLB - CEL*CUSB - P2 (NX) 
R2(3) CRB - (DERBY + P1L(NX)*FVB- (CEL+P2(NX)}*USB + CEL* 

(FV8 + CVB*FB - VB*CFB - CFVB)) 

ENDIF 

R1(J} = F{J-2,2) - F(J,2) » DETA{J-1) *UB 

R3(J-1)= UtJ-1,2) - UlJ,2) * DETA{J-1} *VB 

CONTINUE 


BOUNDARY CONDITIONS 


R1(1} 9.0 
R2 (2) 0.¢ 
R3 (NP) o.c 


RETURN 
END 
COMMON /BLCO/ RL,NBL(2!,XCTRI{2),ntflag,cransnew{2) , NT 


COMMON 
COMMON 
COMMON 
COMMON 
DIMENS 


OPEN 
CPEN 
OPEN 
OPEN 


fBLCL/ 

7BLC2/ 

/BLC3/ 

JBLES/ 
ION 


ITR, XCT 
NX, NXT, 
x(200),0 
DLS {200} 
NXTSF (2) 


xC (200) , YC(200} 
&,NPT,NTR, IT, ISF 
£ {200}, 1 (200), P2 (200), GMTR(200) 
2 V8 (209), CF (200), THT (200) 
2XI (290), ¥1(200}, VEI (200) 


(UNIT=9, FELE 


='bl2d.dat’, STATUS= 


(UNIT=8, FILE='E)2d.out’ , STATUS 


? UNKNORN’ ) 
UNKNOWN’ ) 


(UNIT=20, FILE 
(UNIT=21, FIL 


cf,dat’, STATUS=" UNKNOWN’ } 
=‘ dls.dat’, STATUS=" UNKNOWN’ } 


WRITE(6,*) ‘READING TEE DATA...’ 

READ ( 9,15 ) RL, XCTRI(1), XCTRI(2) 

READ ¢ 9,10 ) NI, IS, ITRANS 

READ ( 9,15 ) (XI(I),¥I{I), VEI(1), T=21,NI) 
WRITE(6,*} ‘INPUT OF DATA COMPLETE.’ 


WRITE(8,90) RL, XCTRI (1), XCTRI(2 
NXTSF(1})= NI - IS + 1 
NXTSF(2)= IS 


DATA FOR EACH SURFACE 
DO 200 ISF = 1,2 

nt £lag=0 

NXT = NXTSF(ISF} 

GO TO (201,202), ISF 


UPPER SURFACE 

II = is-1 

DO 211 I=1,NXT 
II = Irel 
XC(I) = XI(II} 
ye(L) YICII) 
UE{I} = VEI{II) 

CONTINUE 

GO TO 300 





LOWER SURFACE 

Il = Is+1 

DO 212 I=1,NXT 
II = 18-1 
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214 XC(1) = XECIT) 
215 yo(T) = YECRI) 
216 UE(I} = VEI{ID) 
212 CONTINUE 






















xd) = 0.0 
220 DO 301 I=2,NXT 
X(T) = X(T21) +SQRT ( (KC (T) XC (TH 1) 882+ (CCT) “VC (En 2DD 8824 








TRANSITION LOCATION 
225 DO 320 I=1,NXT 

226 GMTR(I) = 0.0 
227 IF (XC{1) .GE, XCTRI{ISF)) GO TO 321 
228 320 CONTINUE 

229 321 NTR = 1 

230 PGAMTR = 1200. 

231 RXNTR = X(NTR-1)**UE(NTR-1) * RL 
RL**2/RXNTR**1. 34*UE (NTR-2) *73 
9.0 

0,5/UE(NTR-1)/ PGAMTR 

235 DO 322 I = NTR, NXT 

0.5/UE (1) /2GAMTR 

UEINTG+ (UL+UZ} * (X(T) HX LEAD) 
U2 
GGFT*UEINTG* (X(1) -X¢NTR-1)) 
10.0} GO TO 323 

241 GMTR(I) =_ 1.0-EXP {-GG) 

242 322 CONTINUE 

243 323 DO 324 IT=3,NXT 

GMTR{II) 1.¢ 









PRESSURE GRADIENT PARAMETERS 
X{2)-XQ) 
UE(2)-VEC) 


ound 








249 ANG2 ATAN2 (DUE, DX? 
250 DL2 DX 

251 DO 331 I = 2,NXT-2 

252 ANGL = ANG2 

253 DLL CL2 

254 DX X{T+1)-%CD) 


255 DUE VECT+1) -UE CL) 


neuen 


256 ANG2 ATAN2 {DUE, DX) 

297 DL2 DX 

258 ANG {DL2*ANG]+DL1°ANGZ} /{DL1+DL2? 
259 P21) TAN (ANG) 

260 331 CONTINUE 

261 P2{NXT) = 2.*DUE/DL2 - P? (NXT-1) 

262 be 330 f 2,NXT 


263 P2(1} = X€X) * P2CT) /VECH) 
264 P1{t) = 0.5 * (1.0 + P2UT}? 
265 330 CONTINUE 
266 P2(1) 
Pitl} 






zl. 
20.5 * (1,0 + F201)} 









270 BOUNDARY LAYER CALCULATION 
27 WRITE(6,*) ‘BOUNDARY LAYER COMPUTATIONS IN PROGRESS...” 

272 CALL BL 

273 WRITE (8, 910) ISF, (EXC (1) X CE), VR ET) CP CD), DLS CT), THT (2), T= 1, NXT} 
274 AECISF.EQ.1) then 

278 write(20,9C5) (XC(I),CF(1},! 
276 write(21,905) (XC€I),DLS(7), 
227 end if 

278 905 FORMAT(FS.4,4x,£11.4) 

CONTINUE 





aNXT} 
NXT) 












281 ***IF AOA is 0 deg., make trans. locs. equal: 
if (vei (2}.eq.vei(ni-1)) transnew(1) =transnew (2) 


{£{ITRANS.eq.0) then 
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print *,’Estimate for upper transition:',transnew(1) 
print *,’Estimate for lower transition:’,transnew(2) 
endif 


CLOSE (UNIT=8} 
CLOSE (UNIT=9} 
STOP 


FORMAT (355) 

FORMAT (3F10.0) 

FORMAT (//5X,‘RL=",£12.5,5X,’XCTRI(1} =',F8.3,5X;"XCTR(2) =", F8.3) 

910 FORMAT(///2X,'*** SUMMARY OF BOUNDARY LAYER SOLUTIONS OF ISF =',12 

+/72X, "NX? , AX, 1XC! BX, 7S", BX, VR? BX, CF’, BX, "DLS! , BX, “THT! 
4/(15,2F8.4,4E11.4)) 

END 

SUBROUTINE EDDY 

COMMON /BLCO/ RL, NBL(2),XCTRI(2), nt flag, transnew(2),NI 

COMMON /BLC2/ NX, NXT,NE, NET, NTR, IT, ISF 

COMMON /BLC3/ X (200) , UE(200) , P1 (200), ®2 (200), GMTR (200) 

COMMON /BLC7/ ETA(201),DETA(201) ,A(201) 

COMMON /BLCB/ F(201,2),U(201,2}, ¥(201, 2), B{201, 2) 

DIMENSICN EDVI (201) 


RL2 SQRT(RL*UE {NX} *X (NX}} 
RLA SQRT(RL2) 
RL216 0.16 * RL2 


ALFA 0.0168 
EDVO ALFA‘RL2*GMTR (NX) *{U(NP, 2) *ETA(NE) -F (NP, 21) 
EDVI{1)= 0.0 
YBAJ RL4* SORT (ABS (V (1,211) /26.0 
DO 70 J=2,NP 

Js = 
YBA BAJ*ETA (J) 

10 


EL 
IF(YBA .LT, 10,0) EL = 1,0 - EXP(-YBA} 
EOVI(J}) = RLZ16*GMTRINX) *(EL*ETA(J}) **2 * ABS(VI5,2)) 
LF{EDVI(J) .GT. EDVO} GO TO 30 
IF (EDVIGS) LE. EDVI(J-1)) EDVI(G)= EDVI(J-1}) 
B(J,2) = 1.0 + EDVI(J) 

CONTINUE 


B(1,2) = 1.0 


RETURN 

END 

SUBROUTINE OUTPUT 

COMMON /BLCO/ RL,NBL(2),XCTRI(2),ntflag,transnew(?!,N1 
COMMON /BLC2/ NX,NXT,NP,NPT,NTR, IT, ISF 

COMMON /BLC3/ X(200}, UE{(200}, 11200), P2 (200) ,GMTRI2E3} 
COMMON /BLC7/ ETA(201),08TA(201),A(201) 

COMMON /BLCS/ F(201,2),0(201,2),V(201,2),8(201,2) 
COMMON /BLCS/ DLS(200), VW(200},CF{200), THT (200) 
dimension rdiff(201},rdlow(201) 


TF(NX.EQ.1 } THEN 
DLS(NX)= 0.0 
THT (NX} 

CF (NX) 
Vi (NX) 
rdi ££low=1000 
nstop 
ELSE 





SQRX SORT (UE (NX) *X (NX) *RL} 

CF (NX) 2.0 * Vél,2) * BC1,2) /SORX 
VW (NX) v(1,2) 

DLS (NX)= X(NX)/SQRX * (ETA(NP)-F{NP,2)) 
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- Ul = UQ1,2) * (1,0 -Ut1,2)) 
SUM 9.0 
DO 20 J=2,NP 
U2 U(d,2) * (1.0 -U(0,2)) 
SUM suM + A(J) * (U1 + U2) 
Ul 2 
“ CONTINUE 
THT(NX}= X(NX}/SQRX * SUM 
rex=UE (NX) *X (NX) *RL 
xtheta=UE (NX) *THT (NX) *RL 
rtrans=1,174* (1.0+22400.0/rex) trex**0.46 
rdi ff (nx) =abs (rtheta-rtrans) 
4£ CONX.gt.2) wand. (rdlow(nx-1) eq. rdlow(nx-2)}) then 
if (rdlow{nx-2) .eq.rdlow{nx-3)} nstop=1 
endif 
if (1SF,eq.2) then 
if ((ntflag.eg.1) .and. (nstop.eq.0)) then 
transnew(ISF) =rex/(RL*UE (NX) ) 
nt flag=0 
endif 
endif 
if ((rdiff€(nx) .LT.rdifflow) .and. {nstop.eq.0}} then 
transnew(ISF)=rex/(RL*UE(NX)) 
rdlfflow=rdiff (nx) 
nt flag=1 
endif 
rdlow(nx) =rdifflow 
ENDIF 
print *,nstop,nx, rdiff{ nx-1),rdiff(nx-2} 
print *, isf,rex,transnew(isf) 
print *,rtheta, rtrans, rdiffinx),rdifflow 


SHIFT PROFILES FOR THE NEXT STATION 


ymark=.0005 
DO 175 J=1,NPT 
1f(ISF.EQ.1) then 
L£(UC2, 1) LT, (0.995) chen 


TA (3) *SQRT(X (NX) / (RL*UE (NX) }) 
do nxloop=5,NI/2-1,5 
if (NX.EQ.nxlocp) then 
markKx=NX/5 
nu ma rkx- 30 
*) OT, 1) emarkx, yplot 
write (60,°) U{J,1)+markx, yplot 
if(yplot.gt.ymark) then 
write (S5,*) markx, ymark 
ydiff=yplot-yplotold 
udiff=U(J,1)-U(J-1,1) 
xvalvesU(J~1, 1) +udiff* (ymark-yplotoid) /ydif 
write(5S,*) xvalue+markx, ymark 
write (55,92) 
format {/} 
ymark=ymark+.0005 
if (yplot.GT.ymark) goto 91 
endif 
endif 
end do 
else 
if (lasty.EQ.1) then 
lasty=0 
do m=1,2 
do nxloop=$,NI/2-1,5 
if (RX.EQ.nxloop) then 
markx=NXx/5 
numw=markx+30 
write (numw,*) markx, yplot 
write (60,*} markx, yplot 
endif 








HS 
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end do 












428 ypiot=0.0 
429 end do 

430 endif 

431 endif 

432 endif 

433 yplovoid-yplot 

434 F(S,1) = F(J,2) 
435 u(s,1) = 0,2) 
436 vQJ,1) = V{5,2) 
437 Bi, 1) = B(5,2) 


CONTINUE 








RETURN 









441 END 

442 SUBROUTINE SOLV3 

443 COMMON /BLC2/ NX,NXT,NP,NPT,NTR, IT, ISF 

444 COMMON /BLC7/ ETA(201),DETA(201},A{201) 

445 COMMON /BLC8/ F(201,2),U{201,2},V(201,2),B(261,2) 

446 COMMON /BLC9/ S1(201),S2 (201) ,53(201),$4{201), $5 (201) ,56(201), 
44? + $7{201) ,$8 (202) ,R1{201) , R2(201) ,R3 (201), R4 (201) 
448 COMMON /BLC6/ DELF{201),DELU(201) , DELV({201} 

449 DIMENSION Al1 (201) ,A12 (201) ,A13 (202) ,A14{201), 


* A21 (201) ,A22 (201) ,A23 (201) ,A24{201) 


464 G22 = G23+S6(2} 

465 AI1(2}= 1.0 

466 A12i2)=-A(2)-G13 

467 A13(2)= A(2}*G13 

468 AZ1(2)= $3(2) 

469 A22(2)= $5(2)-G23 

470 A23(2}= S1{2}+A (2) *G23 

471 R1(2) = R1(2}-(G11*RE {1} +G12*R2 (1) +G13°R3 (10) 
472 R2(2) = RZ(2)-(G21*R1 (1) +G22*R2 (1} +G23*R3(T}} 










c 
Cc 





FORWARD SWEEF 






DO 500 J=2,NP 

477 DEN = (AL3(J=1) *A22 (3-1) -A23 (3-1) *A22 (9-1) -A)” 

478 + {A12 (S-1) 4021 (J-1) -A22 (J-D) AAT (S10 

479 DEN! A22(J-1) *A(J}-A23 (9-1) 

{A23 (J-1} +A(S) * (ACJ) PAZ 1 (5-2) -A22 (2 
¢ 


2 /DEN 
481 G12 = = (AC) *A() 4611? (AL2 (I-12) *A() ~AT3 


ty}} 
139) /DENL 


2 
(GL E*A13 (J-1) +G12*A23 (0-13) /A(S) 

($2 (J) *A2L (T=21-$4 (1 4A23 (I-12) 4A (0) * (54 (0) * 
484 + A22 (J-1} -S6 (J) *A21 (J-1))) /DEN 


4as G22 (-$2 (3) $645) *A (J) -G21* (ACS) 212 (9-2) APB (S-1)}) /DEND 
486 G23 G22*A12 (3-1) +622*A22 (J-1)-SE(T) 

487 All{J 1.0 

48g Al2{s A(J)-G13 

489 Al3(T A(3) *G13 

490 A21(J S$3i3) 

491 A22 (5 $5(3}-G23 

492 A23{T $2{3)+A(I)*G23—- 

493 RJ) R1 (J) - (G12 *R] (J-2) +612 *R2 (I-31) +G13*R3(I-1)) 





494 R2(3) = R2(J)-(G21*R1 (J-1) +G22*R2 (J-1) +G623*R3(J-1)) 
500 CONTINUE 


c 
C BACKWARD SWEEP 


CDELU (NP) R3 (NP) 
El R1{NP)-A12 (NP) *DELU (NP) 
£2 R2 (NP) -A22 (NP) *DELU (NP) 
DELV (NP) (B2*A11 (NP) -EL*A21 (NP) ) /(A23 (NP) *AL1 (NP) ~A13 (NE) * 
+ A21{(NP)) 
DELF (NP} (E1-A13 (NP) *DELV (NP) } /ALL (NP) 
NP-1,1,-1 

= R3 (3) -DELU(J+1) +A (J+1) *DELV (341) 
A21(I} *AL2 (J) *A (941) -A2Z1 (9) FALS(J}-A(J#1) 4022 (o)* 
Atl (J) +823 (J) #ALL Co) 

DELV(J) = (AL1 (3) * (R2 (3) +3*A22 (J) )-AZ1 {J} *R1 (I) -E3°A21 (I) * 

A12(J)) /DEN2 

DELU(3} e-A(J+1) *DELV (J) -E3 

DELF(J) = {R1 (J) -AL2 (3) *DELU (3) ~A13 (J) *DELV (J) ) /AL1 (Jd 
CONTINUE 





DO 700 J=1,NP 
F(3,2 F (J, 2) +DELF (J) 
v(I,2 US, 2) +DELU (J) 
V(d,2)= Vd, 2} +DELV (J) 

CONTINUE 

gtl,2)= 0.0 





RETURN 
END 
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Output file: bI2d.out 


RL= 0.540005+06 XCTRI{1} = 0.380 XCTRI2) = 0.762 


‘“** SUMMARY OF BOUNDARY LAYER SOLUTIONS OF ISF = 1 


z 
x 


XC s VR cr DLs THT 

9.0024 9.0000 0.1232£+01 0.0000E+00 0.0000E+00 0.00C0E+00 
0.0005 0.0057 0.8456E+00 0.6098E-01 0.1815E-03 0.7901£-04 
0.0005 0.0112 0.1242E+02 0.3873E-01 0.12155-03 0.5674E-04 
©.0024 06.0170 0.1083E+01 0.2228F-01 0.1219£-03 0.5399E-04 
©.0064 0.0236 0.8010£+00 0.1292E-01 0.16566-03 0.6935E-04 
0.0223 0.0314 0.6520£+00 0.8809E-02 0.2164E-03 0.8968E-04 
9.0201 0.0406 0.5141E+00 0.6033E-02 0.2828E-03 0.1131E-03 
0.0298 0.0513 0.46398+00 0.4821£-02 0.3479E-03 0.1377E-03 
©.0433 0.0637 0.4023£+00 0.3749E-02 0.4201E-03 0.1631E-03 
0.0547 0.0777 0.3660E+00 0.3096E-02 0.4966£-03 0.1901£-03 
21 (0,0698 0.0934 0.3224E+00 0.2498E-02 0.5820E-03 0.2185E-03 
22 0.0866 0.1106 0.3097E+00 0.2212E-02 0.6613E-03 0.2467E-63 
13° 0.1059 0.1293 0.2804E+00 0.1860E-02 6.74950E-03 0.2755E-03 
14 0.1251 0.1496 0.2644€+00 0.1637E-02 0.8387£-03 0.3051£-03 
15 0.1466 0.2712 0.2426€+00 0.1411E-02 0.9339£-03 0.3353E-03 
16 0.1695 0.1941 0.2393E+00 0.1311E-02 2.1020£-02 0.3646E-03 
27 0,1937 0.2184 0.2091E+00 0.1085E-02 0.1127£-92 0.3954E-03 
18 0.219% 0.2438 0.1926E+00 0.9513E-03 0.12435-02 0.4278E-03 
19° 9.2456 0.2704 0.1784E+00 0,8401E-03 0.1349E-02 0.4593E-03 
2 0.2731 0.2979 0.16088+00 0.7248E-03 0.14639E-02 0.492 7E-03 
21 0.3015 0.3263 0.1398E+00 0.6046E-03 0.1602E-02 0.5247£-03 
22 0.3307 0.3555 0.1202E+00 0.5001£-03 9.1746E-02 0.558CE-03 
23 0.36CS 0.3853 0.7410E-01 0.2978E-03 0.1965E-02 0.5940F-03 
24 ©.3909 9.4158 0.4473E-01 0.1739E-03 S.2177E-02 0.6315E-03 
25 0.4218 0.4467 0.1177E+00 0.4439E-C3 0.19945-02 0.6692E-03 
26 0.4530 0.4778 0,2823£+00 0.1033E-02 0.1746£-02 0.7035E-03 
27 C.4843 9.5092 0.4308E+00 0.1534E-02 €.)655E-02 0.7520F-03 
28 C.5156 0.5407 0.5688E+00 0.1976E-02 0.1642£-C2 0.8127E-03 
29 0.5469 0.5720 0.7246E+00 0,2460E-02 0.1656£-02 0.8794E-03 
30 0.5782 0.6033 0.8574E+00 0,28452-02 C.17005-02 0.9503E-03 
32 (0.6096 0.6342 0.9414E+0C 9.3062E-02 0.1782E-02 ¢.1031E-02 
32 0.6394 0.6648 0.1017E+01 0.3247£-02 0.18815-G2 0.1118E-02 
33° 0.6693 0.6948 0.1077£+01 0.3379E-02 0.19872-02 0.1206E-02 
34 (0.6984 0.7241 0.1119E+01 0.3454£-02 9.2100£-02 0.1296E-02 
3S 0.7268 C.7526 0.2233E+01 0.3443E-02 0.2229E-02 0.13390E-02 
36 0.7544 0.7803 0.1346E+01 0.3438E-02 0.23264£-02 0.24865-02 
37 0.7808 0.8070 0.1152E+01 0.3411£-02 €.2499E-02 C.1581E-02 
38 ©.8062 0.8325 0.1131E+01 0.3317E-02 0.2655£-02 0.1684E-02 
39° 0.8304 0.8569 0.1107E+01 0.3216E-02 0.2813€-02 0.1787£-02 
40 0.8533 0.8799 0.1089£+01 0.3139£-02 €.29755-02 0.1892£-02 
41 0.8748 0.9016 0.1045£+01 0.29898-92 0.31482-C2 0.1999E-02 
42 0.8949 0.9218 0.9844E+00 0.2804E-02 0.3351£-02 0.2119E-02 
43° C.9133 0.9404 0.9361E+00 0.2€55E-02 0.3557£-02 0.2239E-02 
44 0.9301 0.9573 0.8514E+00 0.2408E-C2 0.3787E-02 0.2366E-02 
45 0.9452 0.9726 0.7083E+00 0.2002E-02 6.4079£-02 6.2513E-02 
46 0.9586 0.9860 0.4504E+00 0.1275E-02 0.4490E-02 0.2688E-02 
47 0.9701 0.9977 0.2796£+00 0.7951E-03 C.4909©-02 0.2877E-02 
48 0.9798 1.0075-0.1193£+01-0.3410E-02 0.7082E-02 0.3153E-02 
49 0.9876 1.0154-0.9638E+03-0.2784E+01 0.1712E+03-0.4403E+07 
SO 0.9936 1.0213 0.2408E+12 0.7077E+99-0.5044F-08-0.3973E+19 
$1 0.9975 1.0253 0.3259E+12 0.9806E+09-0.2188E+09-0. 44945420 
52 0.9995 1.0273 0.940GE+12 0.2907E+10-0.3121£+09-0.5061E+20 


a 
OWOd HM NmBwWYY 











*"** SUMMARY OF BOUNDARY LAYER SOLUTIONS OF ISF = 2 


NX xc s vw cr DLs THT 

0.0024 0.0000 0.1232E+01 0.0000E+00 0.0000£:00 0.0000E+00 
0.0064 0.0056 0.8911E+00 0.4236E-01 0.13162-03 0.5759E-04 
0.0123 0.0244 0.7976E+00 0.2157E-01 0.1814E-03 0.7867E-04 
0.0201 0.0236 0.7397E+00 0.1440E-01 0.2262E-03 0.9751E-04 


owe 








0.6680E+00 0.1026E-01 
0.6251E400 0.7990E-02 
0.5723E+00 0,6280E-02 
0.9379E+00 0,5202E-02 
0.5121E+00 0.44268-02 
0.4870E+00 0.3818E-02 
0.4627E400 0.3321£-02 
0.4419E+00 0.2932E-02 

. 0.42836+00 0.2644E-02 
012014 0.4090E+00 0.2364E-02 
0, 3842E+00 0.2092E-02 
0.3742E+00 0.1928E-02 
0,3582E+00 0.1753E-02 
0.3411E+00 0.1593E-02 
C0.3272E+00 0.1462£-02 
0.3032E+00 0.1301£-02 
O.2810E+00 0.1162E-02 
0.2691E+00 0,1075E-02 
0,2643E+00 0.1021E-02 
0.2476E+00 0.9272E~-03 
0.2102E+00 0.7655E-03 
0.1983E+00 0,7038E-03 
0.1997E+00 0.6912E-03 
0.1641E+00 0.5549E-03 
0.12798+00 0.4237E-03 

16778 O.1173E+00 0,3810E-03 
0,9714E~01 0,3098E-03 
C.3892E-01 0.1220E-03 
0.76048-01 0,2349E-03 
0,5406E~02 0,1646E-04 
18155-0.8287E-01-0.2493E-03 
-9041E-01-0,.2691E-03 
.15626+00-0 . 4602E-03 
12511E+00-0.7339E-03 
1 3360£+00-0.9753E-03 


12989E+07-0.8597E+04 
1 3229E+07-0. 92716404 
1 1935E412-0.5553E+09 
19807-0,8049E+14-0.2316E+12 
1990S 0,1412E+17 0,.4082E+14 
9984 NaN NaN 
0044 NaN NaN 
0083 NaN NaN 
0104 NaN NaN 
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0.2790E-03 
0.3328E-03 
0.3921E-03 
0,4550E-03 
0.S160E-03 
0.5802E-03 
0.6472E-03 
0.7162E-03 
0.7934E-03 
0,8539E-03 
0.9338E-03 
0.1008E-02 
0.1084E-02 
0.1166E-02 
0.1247E-02 
0.1341E-02 
0.1442£-02 
0.1$35E-02 
0.1618£-02 
O.1711E-02 
0,18445-02 
0.1962E-02 
0.20436-02 
0.2182£-02 
0,23668-02 
0,2504E-02 
0.2658£-02 
0.2973E-02 
0.2607E-02 
0,30366-02 
0.33906-02 
0.3384E-02 
0.3787E-02 
0.44685-02 
0.5324E-02 





0..7034E+03- 
-0,1S13E+07- 
0.1354E409- 
0.2118E+10- 
0.307084 14- 
NaN 
NaN 
NaN 
NaN 


0.1184E-03 
0.1401E-03 
0.1629E-03 
0.1874E-03 
0,.2112E-03 
0.2358£-03 
0.2609E-03 
0.2867£-03 
0.3124£-03 
0,3378E-03 
0.3653E-03 
0,3923E-03 
0.4191E-03 
0.4469E-03 
0.4743£-03 
0.5032£-03 
0.5333E-03 
0.5629£-03 
0.5909E-03 
0.6187E-03 
0,6497E-03 
0. 6808E-03 
0.7082E-03 
0.7375E~03 
0. 7697E-03 
0.7996E-03 
0.8285E-03 
0.8592E-03 
0.7663E-03 
0.7969E-03 
0.7639E-03 
0.7997E-93 
0.8365E-03 
0.8788E-03 
0.9352E-03 





1 7438E+01-0.21S0£-01 0.2542E-01-0.54216+00 


0.2021E+09 
0.1238E+16 
OQ. 27465 +22 
~85026+23 
0. 7058E+30 
-NaN 
-NaN 
“Nan 
“NaN 





Gnuplot command file: profile 


set terminal tek40xx “ 
set nogrid 

set. nolabel 

set size 1,1 

set data style lines 

set noxtics 

set ytics 

set title "Velocity Profiles AOA=10 deg." 0,0 

set nokey 

set xlabel “Airfoil Upper Surface Station" 0,0 
set xrange [0 : 10) 





set ylabel *y/c” 0,.5 
set yrange (0 : .010) 
set size .6,.6 

set label 
set label 
set label 
set label 
set label 
set label 
set label 
set label 
set label 
set label 
set label 
set label at 
Plot “prolo","Lin10" 


tation 0 is" at .8,.008 
tagnation" at .8,.007 






center 





Input file: FOROO!. DAT (Stardent) or fort.1 (IRIS) 


3 


NACA 001? AIRFOIL 


aaa 


* ALP PIVOT 

6.000000 0.250000 

NLOWER — NUPPER 

50 50 

x/C 
1.00000 0.98000 0.96000 0.94000 0.92000 0.90000 
9.88000 0.86000 9.84000 0.82000 0.89000 0.78000 
0.76000 0.74000 9.72000 0.70000 0.8000 0.66000 
0.64000 0.62000 0.60000 0.58000 0.56000 0,54000 
9.52000 0.50000 0.48000 9.46000 0. 44000 0.42000 
0140000 0.38000 0.36000 0.34000 0. 32000 0.39000 
9.28000 0.26000 0.24000 9.22000 0.20000 0.18000 
$116000 0.14000 «9.12000 0.10000 9.08000 0.06000 
0.04000 0.02000 0.00000 9.02000 0.04000 0.06000 
6.08000 0.10000 «0.12000 0.14000 0.16000 0.18000 
9120000 0.22000 0.24000 0.26000 0.28000 0.30000 
0.32000 0.34000 0.36000 9.38000 0. 40000 0.42000 
0.44000 0.46000 0.48000 9.50000 9.52000 0.54000 
0156000 0.58090 0.60000 0.62000 0.64000 0.66000 
0168000 0.70000 0.72000 9.74000 0.76000 0.73000 
0780000 0.82000 0.84000 0.86000 0.88000 ©. 90000 
0.92000 0.94000 0.26000 0.98000 1.00000 


r/c 
=0,00126 -0.0C403 -€,00674 -0.00938 -0.01196 -0,01448 
Tolo1e94 0.01935 -0.02170 -0.0239¢ 0.02623 -0.02842 
-0.03056 -0.03264 -0.0346? =0.03664 -0.03856 -0.04042 
0.04222 -0.04396 -0.04563 0.04723 -0,.04878 -9.05026 
-0.05265 0.05294 108415 -0.05530 -C.05634 
-9105803 -0,08868 -C.05923 -0,05366 0.05995 
0.05997 -0.05966 -0.05911 -0.05838 -0.05737 
-0.05444 0.05736 -0.94990 0.04683 -0.04399 
0.03231 -0.02382 0.90006 0.02382 0.03231 
0.04309 0.04683 0.04990 0.05236 0.05444 
0.05737 0.95833 0.05911 0.05966 0.05997 . 
0.05995 0.05966 0.05923 0.05868 0.05803 0.95726 
0.05634 9.05530 0.08415 8.05294 9.05165 €.05026 
0.04878 C.04723 0.04863 0.04396 0.04222 C.04042 
0.03856 6.03664 0.03467 0.03056 9.02842 
0.02623 0.02399 9.02170 0.01694 0.01448 
0.01196 0.00938 0.00674 0,00126 








Input file: incompbl.dat 


IWAKE NXT NY ITREND 
1 16] 37 40 
ite ()) ITR(2) ISWPMX RL XCTR(1} 
o e 1 540000.0 0.30000 
Te 
L 








Output file: incompbl.out 


NACA 0012 AIRFOIL ES 


nea 


INPUT DATA FOR INVISCID-FLOW CALCULATIONS 


ALPI= 6.0000 PIVOT= 0.2506 


NLOWER= 50 NUPPER= Fit) 


COORDINATES OF THE BODY 

x/c 
1.000000 0.980000 0.960000 0.940090 0.920000 9.900000 
0.880000 0.860000 0.840000 0.820000 0.800000 0.780000 
0.760000 0.740000 0.720000 0.700000 0.680000 0.660000 
0.640000 0.620000 0.690000 0.580000 0.560000 0.549900 
9.820000 6.500000 0.480000 9.460000 0.440900 0.420000 
9.400000 0.38000C 0.360900 0.349900 0.329000 0.300000 
0.28000 0.260000 ©.240CC2 9.220000 0.200000 0.130090 
0.160000 0.149000 0.120000 0.100000 0,08c0co 0.060000 
0.940000 0.020000 ©.0C0000 0.020000 0.040000 ¢.c60000 
0.080000 0.100000 ¢.120000 0.140000 0.160000 0.180000 
0.200000 0.220000 0.240000 0.260000 0.280000 9.300006 
0.320000 ©.34000C ¢.360000 0.380000 0.400C9C 6.429000 
0.440000 9.460000 0.480000 0.500000 0.520000 9.540000 
0.560900 0.580002 ¢.600000 0.620009 0.640900 9.660000 
O.68N00D H.FCCODD 0.720000 $.7400CO 9.760000 0.780000 
0.800000 0.829C09 C.840C00 O.86C900 C¢.880900 2.900900 
©.920906 0.940000 9.969000 O.9800c0 7.000900 


+i 


=0.001260 -0.004030 -0.006740 -0.009380 -C.011960 -0.014480 
0.016940 -0.019350 -0.621700 -0.023990 -0.026230 -0. 28420 
-0.030560 -0.032640 -5.034670 -0.936640 -0,038560 -0.C40426 
0.942220 -0.043960 -0.045630 -0.047230 -0.048789 -G.0S0260 
-0.051650 -0.052940 -0.054150 -C.055300 -0.05634o -0.957260 
~0.058030 -0.058680 -0.059230 -0.059660 -C.959950 -9.c6006C 
0.059970 -0.059660 -0.089110 -0.058380 -0.95737C -0.05607C 
0.054440 -0.052360 -0.049900 -0.046830 -0.043090 C.038420 
=0.032310 -0.023820 0.0C00CC 0.023820 0.03231¢ 0.038420 
0.043090 0.04683C 0.049900 0,052360 0.054440 0.056070 
0.057370 0.058380 9.959110 0.059660 0.059970 2.060060 
0.059350 0.059660 0.959230 0.058680 0.058030 2057260 
0.056340 0.055300 0.054150 9.052940 0.05265¢ 2.250260 
0.048780 0.04723€ 0.045630 0.043969 0.042220 9.045420 
0.038560 0.036640 0.034670 0.032640 2.930560 ©.928426 
0.026230 0.02399G 0.021700 0.019359 0.526940 6.014480 
0.031960 0.009380 ¢.006740 0.004030 0.001269 
1 INVISCID FLOW RESULTS 





PANEL XP Ye ce 

0.93000E+00 -0.149SCE-02 0.25966E+90 
0.97000E+00 -0.44833E-02 C.18882E+C0 
0.95000£+00 -0.74415£-02 0.14703E+0C 
O.93000E+09 -@.10321E-02 Q.11735E+0C 
0.91000E+00 -0.13061£-01 Q.98371E-01 
2. 89000E+00 -0.15646£-01 0.87319£-01 
0.87000E+00 -0.18113E-01 0.79369E-01 
0.85000E+00 -0.20502E-01 0.72019E-01 
0.83000E+00 -0.22829E-01 0.65720E-02 
0.81000£+00 -0.25104E-01 0,58929E-01 
O.79000E+00 ~0.27325E-01 0.53452E-01 
0.77000E+00 -0.29490E-02 0.48109E-01 
O.75000E+CO -0.31600E-01 0.43794E-01 


URE OCOVHURWNE 








0.73000E+00 
0.71000E+00 
0, 690CCE+00 
0.67009E+00 
0. 65000E+00 
0.63000£+00 
G,61000E+00 
0,59000£+00 
0.57000E+00 
0.55000E+00 
0.53000E+00 
0.51000£+00 
0. 49000E+00 
0.47000E+00 
0.450008+00 
0.43000E+00 
0,41000E+00 
0.39000E+00 
0.37000E+00 
0. 35000E+00 
0. 33000£+00 
0. 31000E+00 
0.29000E+00 
0.27000E+90 
0.25000E+00 
0.23000E+00 
0.21000E+00 
0.19000E+00 
0.17000E+00 
9. 15000E+00 
0.13000E+00 
0,110008+00 
0.90000£-01 
0. 70000E-01 
0. S0000E-01 
0.30000E-01 
0. 10000E-02 
@.100009E-01 
0.30000E-01 
0.50000E-01 
0.70000E-01 
0.9000CE-01 
Q.11L0G0E+00 
0.13000E+00 
O.150C0E+00 
0.17000E+00 
0.19000E+00 
O.21000E+00 
0.23000E+00 
0. 25000E+00 
0.27000E+00 
0,.23000E+00 
0.31000B+00 
0.33000E+00 
0.35000E+00 
0.37000E+00 
0.39000E+00 
0.41000£+00 
0.43000E+00 
0.45000£+00 
0.42000E+00 
©. 49000E+00 
0.51000£+00 
0.53000E+00 
0.S5000E+00 
0.57000E+00 
0.59000E+00 
0.61000£+00 
0. 63000E+00 
0.65000E+00 
0,67000E+00 


70.33655E-03 
-35655E-01 
+37600E-01 
«39490E-01 
-0,42320E-01 
-0.43090E-02 
-0,$4795E-01 
-0.46430E-01 
-0.48005E-01 
-0,49520E-01 
-0.50955E-01 
-0.52295E-01 
-0.53545E-02 
-0,.54725E-01 
-0,55820E~01 
-0.56800E-01 
+S7645E-01 
~0.58355E-01 
-0.589SSE-01 
-0,59445E-01 
-0,5980SE-02 
-0,6000SE-01 
-0.60015E-91 
-0,$9815E-01 
-0,$9385E-01 
-0.58745E-01 
-0,.57875E-01 
-0.56720E-01 
-0.5525SE-01 
~0.53400E-01 
~0.51130E-01 
-0.48365E-01 
-0.44960E-02 
~0.407558-01 
-0.35365E-01 
-0,28065E-01 
-0.11910£-01 
O,11910E-01 
0.28C65E-01 
0.35363£-01 
0.40755E-01 
0.44960E-01 
0.4836SE-01 
0.51130E-01 
0.$34C0E-9) 
0,55255£-01 
0.56720E-01 
O.5787SE-01 
0.58745E-01 
0.59385E-01 
0.59815E-01 
0.6C01SE-01 
0.6000SE-01 
O.598C5E-01 
0.594458-01 
0,58955E-01 
0,5835S5e-01 
0.57645E-C1 
9.56800E-01 
0.55820E-01 
0.5472SE-01 
0.53545£-01 
0,52295£-01 
0.50955E-02 
0.49520E-01 
0. 4800SE-01 
0.46430E-01 
0.44735E-01 
0.43090E-01 
0.413206-01 
0.39490E-01 
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C.39437E-01 
O,35961£-01 
0.32250E-01 
0.29030E~-01 
0.26174E~-01 
0,23331E-01 
0,21248E-01 
0.20161E-01 
0.17719E-01 
O,15010E-0f 
0.13363E-01 
0.13779E-01 
0.14543£-01 
0.12579E-01 
0.11631E-01 
Q,11746E-01 
0,15506£-01 
0.20332E-01 
9,23743E-01 
0,26773E-02 
9.29743E-01 
O.34896E-01 
0.425718-01 
9.53307E-01 
9. 69110E-01 
0.82783E-01 
90.99470E-01 
9.120386+00 
0.14436£+00 
0.18096E+00 
O.21706E+00 
0.26992E+00 
0,33458E+00 
0.42103E+00 
0,54734E+00 
0.70537E+0C 
0.92260E+00 
-0,22412E+01 
-0,228C9E+01 
-0,18259E401 
-0.15937E+01 
-9.14321E+01 
O.13171E+01 
-0,12149E+01 
-0.21409E+01 
0. 10709E+01 
~0,10078E+01 
=0.95329E+00 
~0,89952E+00 
~0.85876E+00 
~0,91972E+00 
~0.78028E+00 
~0.74122E+60 
~0.70214E+00 
0. 66601E+00 
~0,.63339E+00 
0.606145 +00 
-0.58144E+00 
~0.55344E+00 
-0,52502£+00 
-6.49560E+00 
-0,47125£+00 
-0,.44926E+00 
-0.42597E+00 
-0,40167E+00 
-0.37797E+00 
-0,35733E+00 
~0.33679E+00 
-O0.31577£+00 
-0,.29529E+00 
-0,27507E +00 











85 0,69000E+00 0.37600E-01 -0.25461£+00 
86 @.71000£+00 0.35655E-01 -0.23509£+00 * 
87 0.73000E+00 0,33655E-01 -0.21494E+00 
88 9. 75000£+00 0.31600E-01 70,19540£+00 
89 0.77000E+00 0.29490E-01 -0.17526E+06 
90 0.79000E+00 0.27319E-01 -0.15370E+00 
92 0.81000E+00 0.25096E-01 -0.13335E+00 
92 0.83000£+00 0.22823E-01 = -0..11255E+00 ™ 
93 0.85000E+00 0.20492E-01 -0.91752E-0) 
94 0.87000E+00 Q.18084E-01 -0.69232E-01 
95 0.89000E+00 0.15564E-01 ~0.40632E-01 
96 0.91000E+00 O.12899E-01 9 =0.29571E-02 
97 0.93000E+00 0.10105€-01 0. 42584E-01 
98 0,95000E+00 0.72366E-02 0.94521E-01 
99 0.97000E+00 0.43440E-02 0.15781£+00 
100 0,99000£+00 0,14480E-02 0,25966E+00 
1 INVISCID WAKE RESULTS 
PANEL XP YP ce 
101 0.10069E+01 0.839528-04 0.302005+00 
102 0,10222£+01 €.30601£-03 0.21875£+0C 
103 9.104132+01 0. 65883E-03 0.17465E+00 
104 0.10649E+01 0.119438-02 0.14324E+00 
105 0.10942E-01 O.19B21E-02 0.21851£+00 
106 0.113CSE+01 0.31169E-02 0.98103E-02 
107 O.117552+01 0.47269E-92 0.80872E-01 
108 0.123132+01 0.69828E-02 0.66183E-01 
1c9 O.13005E+01 0.10111E-01 0.53640£-01 
110 O.138636+01 0.14408E-01 0,42979£-01 
111 0.14926E+01 0.20260E-01 0,33995E-01 
112 0.16244E+03 0.28163E-01 0.26S15E-01 
113 9.17878E+01 0.38747E-01 0.20374E-01 
114 O.199C3E+01 ©.528448-C1 0.15419e-02 
115 0.22414E:0) «T136SE-01 O.11487£-01 
116 0.255265+01 -956S3E-01 0.84237E£-02 
1:7 0,.29383E+01 +12723£+00 0.6C806E-C2 
118 O.34165E+01 + 16802E+00 0.43224E-02 
119 0. 383646601 «204802400 0.33384E~02 


929000 


CL = 0.71990E+00 


INPUT DATA FOR SOUNDARY-LAYER CALCULATIONS 


IWAKE NXT NW ITREND 
1 16} 37 4c 
ITR() ITR{2) ISWPMX 10**-6*RL XCTROL?) 
0 ie) 1 0.54 0.01 
ie 


1 
OITERATIONS EXCEEDED MAX IN WAKE 8.L. CALCULATIONS AT 
OINTERA= ANX= 1621T= 20 
OCALCULATIONS ARE ALLOWED TO CONTINUE 


Ananda ee hAeee CYCLE dQ tetetneneannee 


BOUNDARY LAYER PROPERTIES FOR THE LAST CYCLE 
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w--+- UPPER SURFACE ---~"~ 


XCTR= 0.850E-01 


NX x 


14 9,00888 
75 0.00540 
76 9,00294 
27° -9,00140 
78 0,00056 
73° 0.00018 
80 0.00003 
81 9.00000 
82 0.00003 
83 0,00018 
84 0.00056 
85 9.00140 
86 6,00294 
87 0.00540 
88 0.00888 
89 9.01333 
990 66.0184] 
91 0.02387 
92 0.03052 
93 0.03806 
94 9.04609 
95 0.05484 
$6 0.06424 
97 0.07430 
98 9,08498 
99 0.09627 
10 9.10814 
101 0.12058 
102 0.33357 
303 0.14708 
104 0.16108 
105 0.17557 
196 0.19052 
107 0,20589 
108 0,22168 
109 0.23786 
110 0.25439 
1110-27127 
112 9.28846 
113. 0.30593 
114° 0.32367 
11S 0.34164 
116 0.35980 
117 (0.37814 
318 0.39667 
119° 0.41535 
120 0.43414 
121 0.45303 
122 0.497197 
123 0.49095 
A124 9.50994 
125 0.52891 
126 0.54783 
127 0.96668 
128 0.58542 
129 0.60405 
130 0.62252 
131 0.64081 
132 0.65890 
133 0.67677 


xs 


9.003702 
9.008237 
0.012084 
0.015235 
0.037694 
0.019500 
0.020776 
0.021748 
0.022721 
0.023996 
0.025802 
0.028261 
0.031413 
0.035260 
0.039795 
0.045011 
0.050900 
0.057456 
0.064668 
0.072526 
0.081020 
©.090139 
9.099869 
0.110198 
0.121120 
0,132591 
0.144626 
0.157198 
0.170288 
0.183882 
0.197955 
0.212494 
0.227475 
0.242880 
0.258685 
0.274871 
Q.291414 
0.308292 
0.325482 
0.342959 
0.360639 
0.378667 
0.396830 
0.415179 
0.433716 
0.452411 
0.471223 
0.490231 
6.509108 
0.528118 
9.547147 
0.566358 
0.585126 
0.604025 
0.622828 
0.642511 
0.660047 
0.678411 
0.696577 
0.714521 


Cr 


0.15032 
9.08157 
0.04906 
0.03732 
0.03117 
9.02762 
0.02532 
0.02382 
© .02238 
0.02073 
0.01865 
0.01614 
0.01327 
0.01023 
0.00757 
0,00579 
0.00474 
0.003294 
0.00300 
0.00184 
0.00082 
0.00026 
9.00004 
0,00026 
0.00137 
0.00306 
0.00438 
9.00493 
2.00507 
9.00497 
0.00481 
0.00466 
0.00453 
0.00442 
0.00431 
0.00421 
Q.9c413 
0.00404 
0.00334 
0.00384 
0.00375 
0.0036? 
0.00361 
9.00356 
0.00351 
0.00345 
0.00338 
0.00330 
0.00324 
9,00320 
9.00316 
0.00310 
0.00304 
0.00298 
0.00294 
&.00289 
0.00285 
0.00279 
0.00274 
0,00269 


DLS 


0.00014 
9.00013 
0.00012 
9.00912 
Q,00011 
0.00011 
0.c001) 
0.00011 
9.00011 
0,00011 
0.00011 
9.00012 
0.00013 
0,00915 
0.90019 
0.00023 
0.00026 
0.00030 
9.00035 
0.00044 
9.00059 
Q.00076 
0.00089 
0.00093 
0.09082 
9.00075 
0.00075 
0.90978 
0.69083 
0.00099 
0.00099 
0.00108 
0.00117 
0.00126 
0.00136 
0.00146 
9.00155 
G.G0165 
0.00176 
0.90187 
0.00198 
0.90210 
9.00222 
0.00234 
0.00245 
@.00256 
0.00269 
0.00282 
0.00296 
0.00309 
0.00322 
0.00335 
0.00359 
0.00364 
0.00379 
0.00393 
0,00408 
0.00423 
0.00438 
0.00454 


125 


UE 


0.13992 
0.32881 
0.52154 
0.70635 
0.86148 
0.97807 
1.06098 
1.12282 
1.18298 
1.25953 
1.36203 
1.48647 
1.61299 
4.70750 
1.75175 
1.76513 
4.76972 
1.77030 
1.76054 
1.73169 
1169096 
1.65398 
1.62448 
1,58915 
1,$3250 
4.49382 
1.47910 
1.43268 
1146562 
1.45392 
1.43861 
1.42270 
1.40737 
1.39279 
1.37837 
1.36513 
1.35299 
4.34127 
1.32926 
1.31691 
1.30442 
4.23247 
1.28062 
1126995 
1.26014 
1.25058 
1.24058 
4.23017 
1.22986 
1.21022 
1.20121 
1.19232 
1.38320 
1.37413 
1.16549 
V.15718 
1.14896 
1.14078 
3.13270 
1117481 


ce iT 


0.98042 
0.89189 
0.72800 
0.50107 
0.25786 
0,04338 
0.12567 
0.26072 
0.39937 
-0.58640 
-0.85513 
-1.20958 
-1.60172 
1.91554 
-2.06861 
2.11862 
=2..43292 
2.13396 
2.09950 
1.99875 
-1.85935 
1.73564 
1.63892 
-1.52839 
1.34548 
-1,23149 
2.18775 
1.16879 
-1.14805 
n1,11389 
1.06960 
1.07407 
-6, 98089 
-9.93987 
0.89991 
0.86358 
~6.83099 
-0.79901 
0.76693 
0.73426 
-0,70151 
066971 
-0.63999 
0.61277 
0.58796 
-0.563294 
0.93903 
-0.91331 
9.48805 
0.46463 
-0.44292 
0.42159 
-0,39997 
~0,37859 
~0.35836 
~0. 33906 
-0,32011 
-0.30139 
-.28302 
-0.26539 


2 
3 
3 
2 
2 
Q 
3 
3 
3 
3 
3 
3 
3 
3 


3 
ES 
3 
3 
3 
4 
4 
5 


4 
9 
2 
4 
3 
2 
2 
2 
2 
2 
2 
2 
2 
2 
1 
2 
2 
2 
2 
2 
2 
1 
1 
2 
2 
2 
2 
1 
1 
2 
2 
2 
2 
2 
2 
- 
2 
2 





we--- LOWER SURFACE 


0.69439 
0.71173 
0.72878 
0.74550 
0.76188 
0.77789 
0.79351 
0.80873 
0.8235) 
0.83785 
0.85371 
0.86503 
0.87795 
0.89029 
0.90208 
0.91331 
0.92397 
0.93404 
0.94351 
0.95237 
0.96061 
0.96821 
0.97517 
0.98147 
0.98712 
0.99209 
0.99639 
1.00000 


0.732218 
0.749643 
0.766773 
0.783584 
©.800052 
0.816156 
0.831873 
0.847182 
0.862061 
0.876489 
0.890448 
0.903918 
0.916880 
0.929317 
0.941212 
0.952546 
0.963307 
0.973478 
0.983046 
0.991998 
1.000321 
1.098003 
1.015035 
1.021406 
1.027108 
2.032134 
1.036475 
1.040126 


XCTR= 0.9446+00 


NX 


xe 


0.01333 
0.01841 
0.02387 
0.03052 
0.03806 
0.04609 
0.05484 
0.06424 
©.07430 
0.08498 
0.09627 
9.10814 
0.12058 
0.13357 
0.14708 
0.16108 
0.17557 
0.29052 
0.20589 
0.22168 
0.23786 
0.25439 
0.27127 
0.28846 
0,30593 
0.32367 
0.34164 
0.35980 
0.37814 
0.39667 
0.41535 
0.43414 
0.45303 
0.47198 


xs 


0.001514 
0.007404 
9.613959 
0.921171 
0.029030 
0.037524 
0.046643 
0.956373 
0.066701 
C.077614 
0.089095 
0.101130 
0.113701 
0.126792 
0.140385 
0.154459 
0.168998 
0.183979 
0.199384 
0.215190 
0.231375 
0.247918 
0.264737 
0.281986 
0.299464 
0.317203 
0.339171 
0.353334 
0.371683 
0.390221 
0.408916 
0.427728 
0.446636 
0.465613 


0.00264 
0.00259 
0.00253 
0.00248 
0.00242 
0.00235 
0.00228 
0.00222 
0.00215 
0.90207 
0.00199 
9.00189 
0.00177 
0.00162 
0.00147 
0.00132 
0.06115 
0.00100 
0.00087 
0.00075 
0.00064 
0.00084 
0.00044 
9.00036 
0.00029 
0.00023 
0.90017 
9.00011 


Cr 


0.35755 
0.13898 
0.04285 
0.02499 
0.01728 
0.01352 
0.01117 
0.00957 
9.00837 
0.00734 
0.00653 
0.00590 
0.00536 
9.00489 
0.00449 
0.00417 
0.00387 
0.00359 
0.00334 
9.00314 
0.00295 
0.00281 
0.00267 
0.00252 
0.00237 
0.00222 
0.00210 
9.00200 
0.00193 
0.00197 
0.00182 
0.00173 
0.00363 
0.00154 


0.00470 
0.00487 
0.00504 
0,00522 
0.00539 
0.00558 
0.00578 
0.00598 
0.00619 
0.00641 
0.00663 
0.00688 
9.00715 
0.00745 
0.00778 
0.00815 
0.00854 
6.00895 
0.90937 
0.00979 
0.01022 
0.01062 
0.01102 
0.01140 
0.01176 
0.01209 
0.03238 
0.01265 


DLS 


0.00014 
0.00016 
0.00038 
0.90021 
0.00026 
0.00030 
0.00034 
0.00037 
0.00040 
0.00044 
0.00048 
0.00051 
0.00054 
0.00052 
0.00062 
0.00065 
0.00069 
0.00073 
0.00077 
0.00082 
0.00084 
0.00088 
0.00092 
0.00095 
0.00100 
0.00104 
0.00209 
0.00113 
0.00117 
0.00121 
0.00124 


0.00328 


0.00133 
0.00138 


126 


1,11693 
1,10914 
1.10245 
1.09380 
1.08618 
1.07852 
1.07085 
1.06328 
1.05581 
1.04838 
1.04096 
1.03340 
3.02559 
1.02739 
1.00885 
1.00016 
0.99160 
0.98336 
0.97557 
0.96835 
0.96173 
0.95571 
0.95029 
0.94548 
0.94126 
0.93761 
9.93453 
0.93195 


UE 


0.05723 
0.26180 
0.44075 
0.$6141 
0, 62934 
0.67535 
0.71733 
0.75644 
0.79113 
0.81986 
0.84452 
0.86609 
0.88542 
0.90187 
0.91609 
0.92896 
0.94062 
0.95018 
0.95838 
0.96566 
0.97198 
0.97777 
0.98337 
0.98811 
0.99165 
0.99417 
0.99595 
0.99732 
0.99871 
1.00031 
1.00197 
1.00311 
1.00340 
1,00306 


-0,24753 
0.23019 
0.21318 
0.19640 
-0,17979 
0.16328 
0.44672 
0.13056 
0.11473 
-0,09910 
~0.08359 
0.06792 
-0.05184 
0.03507 
0.01779 
0.00032 
0.01672 
0.03239 
0.04826 
0.06229 
0.07507 
0.08662 
0.09694 
0.10606 
0.43403 
0.12088 
0.12666 
©.13147 


cP 


0.99672 
0.93146 
0.80574 
0.68482 
0.60393 
0.54390 
9.48543 
0.42780 
0.37414 
€.32783 
0.28679 
0.24990 
0.21604 
0.18664 
0.16095 
9.13704 
9.11524 
0.09716 
0.08151 
0.06751 
0.05526 
0.04396 
0.03298 
0.02363 
0.01662 
0.01163 
0.00808 
0.00535 
9.00258 
-0.90062 
0.00395 
-0.00624 
-0.00682 
0.00612 
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0.49095 
0.90994 
9.92891 
0.54783 
0.96668 
0.58543 
0.60405 
0.62252 
0.64081 
0.65831 
0.67677 
0.69439 
0.71173 
0.72878 
0.74930 
0.76188 
O.77789 
0.79352 
0.80873 
0.82352 
0.83785 
0.85172 
0.86509 
O.87796 
0.89030 
0.90210 
0.91334 
0.92400 
0.93407 
0.94354 
9.9524c 
0.96063 
0.96823 
0.97518 
0.98149 
9.93713 
0.99210 
0.99639 
1.00000 


xs 


1.04119 
1,04769 
1.05492 
1.96297 
1.07193 
1.08191 
1.09303 
1.10541 
1.11919 
1.13454 
.19163 
1.17066 
1.19185 
7.21545 
1.24273 
1.27099 
1.3035? 
4.33985 
1.38025 
1.42523 
1.47533 
1.93212 
1.59322 
1.66239 
1.73941 





0.484624 
9.503652 
0.522663 
0.541632 
0.560530 
0,579334 
0,598017 
0.616553 
0.634916 
0653083 
0.671027 
0.688723 
0.706249 
0.723279 
0.780090 
0.756558 
Q.2772662 
0.788379 
0.903688 
0.818567 
0.832996 
0.846994 
0.860424 
0.873387 
0.885823 
0.897717 
0.909052 
0.919813 
0.929989 
9.939583 
0.948504 
0.986827 
0.964510 
0.971542 
0.977913 
0.983615 
o.988640 
6.992982 
0.996632 


0.92965 
0.91681 
0.92128 
0.92772 
Q.9319¢ 
0.93520 
0.93674 
0.33908 
0.94130 
0.94431 
0.94776 
0.95176 
0.95609 
0.96029 
0.96418 
0.96780 
0.97110 
8.97422 
0.97700 
0.97960 
0.98193 
0.98409 
9.98602 
0.98778 
0.98935 


i) 
o 
q 





5° 
0 
0 
9 
ie} 
Q 
i) 
9 
Q 
0 
0 
Ly} 
Q 
i) 
to) 
0 
c 
0 
9 
Q 
O 
0 
°o 
Qo 
e 
-0 
-0 
0 
-0 
-0 
-0 
“9 
~0 
-0 
-o 


° 
re} 
0 
9 
9 
0 
co) 
6 
is} 
vo) 
° 
0 
0 
oO 
to) 
id 
oO 
Q 
0 
ie) 
0 
to) 
Q 
o 
) 


00149 
-00145 
90140 
,00134 
00128 
00123 
-00120 
00116 
-00112 
300108 
-00104 
~00100 
00096 
00092 
-00088 
.00084 
00079 
00074 
00069 
200064 
00059 
-00053 
~00047 
00041 
~00034 
00026 
00018 
-00009 
-00002 
90006 
- 00016 
-00030 
~00950 
-00074 
209106 
200144 
200181 
00208 
~00223 





ce 


224318 
«15945 
~1S8124 
13935 
-13157 
212541 
-12251 
11813 
~21394 
-10828 
1037S 
-09415 
08589 
07784 
07035 
+0633? 
«05696 
-9§030 
-04947 
-04038 
~93582 
-03157 
202777 
92428 
02118 


9.00143 
0.90147 
g,00150 
0.00155 
9.00260 





0.00523 


DLS 


0.01782 
0.01604 
0.01449 
0.03336 
0.01264 
0.01208 
0.01189 
0.03128 
0.01092 
0.01051 
O.03012 
0.00974 
0.00937 
0.00903 
0.00873 
0, 00846 
0.00822 
g.gteoe 
0.00781 
0.00763 
0.00748 
0.00733 
9.00720 
9.00709 
0.00698 
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1.00262 
2.00243 
1.00230 
1,00178 
1.00087 
0.99985 
0.99895 
9.99806 
0.99702 
0.9958) 
0.99450 
0.99310 
0, 99162 
0.99006 
0.38842 
0.98668 
0.98483 
9.98287 
0.98080 
0.97862 
0.97637 
0.97407 
0.97173 
0.96935 
0.96693 
0.96448 
0.96204 
0.95970 
Q.95757 
0.95569 
0.95402 
0.9823) 
0.95021 
0.99748 
0.94382 
0.93889 
0.93267 
0.92576 
0.93935 


2.33014 
2.00843 
1.93997 
1.79191 
1.68424 
163208 
1.58876 
1.94734 
4.50866 
1.47057 
143422 
1.39982 
1.36780 





0.00924 
=0,00486 
9.00460 
-0.00357 
0.00173 
9.60030 
0.00210 
0.00387 
0.00595 
0.00836 
0.01096 
0.01375 
9.03670 
0,01978 
0.02303 
0.02647 
0.03011 
9.03397 
0,63804 
0.04232 
0.04670 
0.95118 
0.05573 
0.06037 
0.96505 
0.06978 
9.07448 
0.97897 
0.98306 
0.08665 
Q. 08985 
0.09310 
©,09711 
0.40229 
0.10921 
Q.11843 
0.13013 
0.34296 
0.15480 


3 
2 
2 











186 
187 
188 
189 
190 
191 
192 
193 
194 
195 
196 
19? 
198 


1.8251? 0.99078 
1.92067 0.99204 
2.02702 0.99317 
2.14543 0.99416 
2.27730 0.99504 
2.42413 0.99580 
2.58763 0.99647 
2.76969 0.99705 
2.97243 0.99755 
3.19818 0.99797 
3.44956 0.99833 
3.72949 0.99863 
4.04119 0.99859 


SUMMARY OF THE DRAG, 
COEFFICIENTS WITH THE CYCLE 


0.01835 0.0cE8? 
0.01585 0.00681 
0.01361 0.00673 
0.01264 0.90666 
0.00990 6.00660 
0.00839 0.00654 
0.00705 0.00649 
0.00590 0.00645 
0.00490 0.00641 
0.00405 0.00638 
0.00333 0.00635 
©,00273 0.00632 
0.00281 0.00000 


CD IS EVALUATED FROM FAR-WAKE FORMULA 
CL & CM FROM INTEGRATION OF CP 


cre! 


VRUSUN EH OO OIRVSWN = OOBURUBLNH AD 


RMR RRNA ee 


cD 
0.012592 


cL cM 


011899 0.665033 
.012128 0.665977 
-011455 9.665906 
-012898 0.665108 
-012117 0.666036 
+011452 0.666010 
011898 0.665143 


0.012272 0.912116 0.666060 


COMPTON DCDODDOCOCCCCOODOSO 


011606 0.011458 0.666076 
012963 0.021896 0.665121 
012167 0.659153 6.002399 
011474 0.658606 0.002286 
011860 0.658523 9.902802 
022045 0.658816 0.002517 
011458 0.657755 ¢.0025¢0 
011853 0.658142 0.002906 
012176 0.658669 9:0025S0 
011600 0.659903 0.002041 
911900 0.661304 0.002273 
012125 0.662937 0.001680 
012497 0.663051 0.001382 
011887 0.663583 0.001815 
9 
oO 
Q 
° 
0 
9 
0. 
ts) 


79 0.001284 


11466 0.664628 0.001075 
11872 0.664564 0.001612 
12094 0.665672 0.001124 
11457 0.665504 0.000910 
11898 0.664903 0.00152¢ 
12119 0.665869 0.001081 
11455 0.665745 0.000863 
011899 0.665033 9.001489 


0.012218 0.665977 0.00106) 
0.011455 0.665906 0.000835 
0.011898 0.665108 0.001471 
0.012117 0.666036 0©.0C1050 
0.011452 0.666010 0.000817 
0.011898 0.665143 0.001461 
0.012116 0.666060 0.001046 
0.011458 0.666076 C.0008C6 
0.011896 0.665151 0.00145? 
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1.12864 
1,1291C 
1.21030 
1.10217 
1.09469 
1.08789 
1.08172 
2.07607 
1.07090 
1.06616 
2.06177 
1.05778 
0.00000 


LIFT AND PITCHING MOMENT 


ONAN NANA NNNN 


Source code: pres.f 


Page 1 


program pres 
REAL cp(100} 
OPEN(unit=20, file='cp.in’, status=! unknown! } 


print *,’no. points before bubble?’ 
read *,N 





print *,/no. points after bubble?’ 
read *,M 





ep(1) corresponds to stagnation for M=.3, dark fringe 5.5 
cp(2) corresponds to next dark fringe 4.5, etc. 


do i=t,N 
read(20,*) cp(i) 
end do 


do i=1,N 
read(1,*) k,x.y 
write(2,*) x,-cp({) 
end do 


do L=N,N-M+1,-1 
read(1,*) k,x,y 
write (2,*°) x,-ee(i) 
end do 


end 
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APPENDIX B: BASIC COMPUTER COMMANDS 


The following collection of commands is intended as a very basic user’s guide to 
the various support programs which are necessary to do research on a UNIX system. 
It is by no means a comprehensive list, simply enough to get started without wasting a 
great deal of time on finding elementary procedures and syntax. In some commands, the 
arbitrary word filename or abbreviation fn is used. In others, specific examples are used 
when it aids the clarity of the explanation. Both methods, however, indicate that the user 


may substitute an appropriate name. 
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man filename 


cp fnl fn2 


ep /alpha/nowak/bl/fn1 . 


cp ../fnl . 
mv fni fn2 
more fn 

Iv 

cd paneldir 
rm fn 


mkdir paneldir 
rmdir paneldir 


Is 
pwd 


batch < fn 


BASIC UNIX COMMANDS 


obtain on-line help information for a program or command 
from the on-line manual 


copy filename! to filename2, both in current directory 


copy fni from another directory, path specified, to the = 
name in the current directory 


copy fni from the directory above to the same name 
move, or rename, fnl to fn2; fn1 will no longer exist 
type the text of the file on the screen, read only 

repeat the last command that started with v 

change directory to paneldir 

delete (remove) fn 


make directory (example name paneldir) 
remove (delete) directory 


list contents of directory (like dir on a pc) 
print working directory 


execute a command file in batch mode (runs even after 
logging off) - useful for long run-time programs 


telnet 131.120.254.92 (suzqt -Stardent) 
131.120.254.91 (madmax - IRIS) 
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BASIC EDITING COMMANDS FOR THE VI EDITOR 


NOTE: ALL commands in VI are case sensitive, type exactly as shown. Check the 
status of the CAPS LOCK key if a command does not seem to work properly. 


vi fal (invokes editor, calls old file if it exists, otherwise creates new 
file) 


To get started into text mode from command mode: 


a (add to document, cursor moves to right and input is enabled) 
i (insert, cursor does not move and input is enabled) 
fe) (open a new line below the cursor, input is enabled) 


To return to command mode: 


Esc (disables input, enables move, write, save, etc.) 


NOTE: All of the following commands assume command mode 


To move around: 


ctrl-h move | left 
ctrl-1 move | right 
ctrl-k move | up 
ctrl-1 move | down 


The above commands will always work. On some machines, the arrow keys MAY work 
as well. Other commands: 


1(shift})G go to first line 

(shift)G go to last line 

/o12D search for the next text string “bl2D" after the cursor location, 
CASE SENSITIVE 


/Cnter key) search for another occurrence of the previous search string ¢ 





When done editing or viewing: 


(Esc): w write, or save, but do not exit 
i (Esc):q quit, exits only if no modifications were made 

(Esc):q! quit, exits without saving even if modifications were made 

(Esc): wq write quit, exits and saves all modifications to the original fn 
. (Esc): wq newfn write quit, exits and saves all modifications to a new fn 

To delete: 

dd delete current line 

2dd delete current line and next line 

10dd delete current line and next 9 lines 


To cut (copy) and paste: 


5yy yanks 5 lines to buffer (leaves original 5 lines also) 
move cursor to desired location 
p pastes the 5 lines 
2 
’ 
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BASIC GNUPLOT COMMANDS 


NOTE: The following is a list of some example commands. Extensive on-line help is 
available by typing help while in gnuplot. 


set term tek40xx (Stardent) 
plot "VEL.DAT" 

plot "VEL.DAT” with lines 

set xrange [0:1] 

set yrange [0:5] 


plot “cf.dat", "dis.dat" (two different data files) 
plot “cf.dat" with lines, "dls.dat" with lines 
plot "bl.dat” using 2:5 (one data file with multiple columns) 


pilot "bi.dat” using 2:5, “bl.dat" using 2:6 

set data style lines (option with lines will then not be needed after each plot command) 
set key (legend) 

set nokey (no legend) 

set grid 

set nogrid 

set nozeroaxis 


To print (these commands can be put in a command file): 
set term postscript 

set output “gnuout" 

replot . 

set term tek40xx (reset to terminal being used) 
set output 
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ftp suzqt 
ftp madmax 


get fn 

get oldfn newfn 
put fn 

quit 


FILE TRANSFER USING FTP 


(numerical computer address may be used instead) 
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